ABSTRACT The conventional 3D magnetotelluric (MT) forward modeling and inversions generally assume an isotropic earth model. However, wrong results can be obtained when using an isotropic model to interpret the data influenced by the anisotropy. To effectively model and recover the earth structures including anisotropy, we develop a 3D MT inversion framework for a triaxial anisotropic model. We use the unstructured finite-element method for our forward modeling. This offers more possibility to simulate more complex underground geology and topography. To solve the inverse modeling problem, we use a limited-memory quasi-Newton algorithm (L-BFGS) with a parallel direct solver for optimization that avoids the explicit calculation of the Hessian matrix and saves the memory and computational time. We validate our algorithm via numerical experiments on both synthetic data and MT survey data from the US Array project.