Three-dimensional (3D) geophysical inverse problems by their very nature are ill posed because of insufficient data and their non-linear characteristics. They also tend to be large in scale. The ill posed nature of these problems can be reduced using regularization techniques with constraints, but the large scale difficulty remains; thousands to tens of thousands of model parameters must be estimated in order to render an accurate map of the subsurface physical properties. In such inverse problems it can be computationally prohibitive to form and invert, repeatedly, a full least-squares system matrix to determine the model update. To avoid these difficulties, iterative gradient methods can be employed. In this paper we will explore the non-linear conjugate gradient (NLCG) method to solve large scale geophysical inverse problems. The application will be directed at determining the subsurface 3D electrical conductivity structure by inverting measurements of magnetotelluric (MT) fields at the earth s surface. These fields arise from the interaction of the solar wind with the earth s magnetosphere. The inverse problem is set up by first forming an objective functional, which is a weighted squared error measure between the observed and predicted data. In order to minimize the objective functional, its gradient with respect to the model parameters is required at each relaxation step of the NLCG scheme. A univariate line search is also required to determine the optimum model update. Preconditioning based on diagonal scaling is further employed to reduce the number of relaxation steps. Here we examine the use of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) secant update to iteratively recur the diagonal of the Hessian matrix of the objective functional. Reduction in the number of relaxation steps with preconditioning can be significant; time savings of a factor of two to three has been observed. In order to compute the predicted data and gradients, finite difference methods have been employed to model the 3D MT fields. Finite difference methods are ideal for computing the MT response arising from complex 3D geological media due to their efficiency. To solve for the fields, a sparse linear system of upwards of a million field unknowns is constructed. Once again because of the size of the system, direct matrix inversion is prohibitive. Instead we employ iterative Krylov solvers to obtain the solution to a predetermined error level. Because the linear system is complex symmetric we use a preconditioned quasi-minimum residual (QMR) solver for efficiency and stability. Simple Jacobi scaling is shown to be an effective preconditioner. * Sandia National Laboratories is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under Contract DE-AC04-94AL85000.