Randall L. Mackie

Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion

We investigate a new algorithm for computing regularized solutions of the 2-D magnetotelluric inverse problem. The algorithm employs a nonlinear conjugate gradients (NLCG) scheme to minimize an objective function that penalizes data residuals and second spatial derivatives of resistivity. We compare this algorithm theoretically and numerically to two previous algorithms for constructing such “minimum‐structure” models: the Gauss‐Newton method, which solves a sequence of linearized inverse problems and has been the standard approach to nonlinear inversion in geophysics, and an algorithm due to Mackie and Madden, which solves a sequence of linearized inverse problems incompletely using a (linear) conjugate gradients technique. Numerical experiments involving synthetic and field data indicate that the two algorithms based on conjugate gradients (NLCG and Mackie‐Madden) are more efficient than the Gauss‐Newton algorithm in terms of both computer memory requirements and CPU time needed to find accurate solutio...

Three‐dimensional electromagnetic modeling using finite difference equations: The magnetotelluric example

We have developed a robust and efficient finite difference algorithm for computing the magnetotelluric response of general three-dimensional (3-D) models using the minimum residual relaxation method. The difference equations that we solve are second order in H and are derived from the integral forms of Maxwells equations on a staggered grid. The boundary H field values are obtained from two-dimensional transverse magnetic mode calculations for the vertical planes in the 3-D model. An incomplete Cholesky decomposition of the diagonal subblocks of the coefficient matrix is used as a preconditioner, and corrections are made to the H fields every few iterations to ensure there are no H divergences in the solution. For a plane wave source field, this algorithm reduces the errors in the H field for simple 3-D models to around the 0.01% level compared to their fully converged values in a modest number of iterations, taking only a few minutes of computation time on our desktop workstation. The E fields can then be determined from discretized versions of the curl of H equations.

Three-dimensional magnetotelluric modeling using difference equations­ Theory and comparisons to integral equation solutions

We have developed an algorithm for computing the magnetotelluric response of three‐dimensional (3-D) earth models. It is a difference equation algorithm that is based on the integral forms of Maxwell’s equations rather than the differential forms. This formulation does not require approximating derivatives of earth properties or electromagnetic fields, as happens when using the second‐order vector diffusion equation. Rather, one must determine how averages are to be computed. Side boundary values for the H fields are obtained from putting two‐dimensional (2-D) slices of the model into a larger‐scale 2-D model and solving for the fields at the 3-D boundary positions. To solve the 3-D system of equations, we propagate an impedance matrix, which relates all the horizontal E fields in a layer to all the horizontal H fields in that same layer, up through the earth model. Applying a plane‐wave source condition and the side boundary H field values allows us to solve for the unknown fields within the model. The r...

Crustal reworking at Nanga Parbat, Pakistan: Metamorphic consequences of thermal‐mechanical coupling facilitated by erosion

Within the syntaxial bends of the India-Asia collision the Himalaya terminate abruptly in a pair of metamorphic massifs. Nanga Parbat in the west and Namche Barwa in the east are actively deforming antiformal domes which expose Quaternary metamorphic rocks and granites. The massifs are transected by major Himalayan rivers (Indus and Tsangpo) and are loci of deep and rapid exhumation. On the basis of velocity and attenuation tomography and microseismic, magnetotelluric, geochronological, petrological, structural, and geomorphic data we have collected at Nanga Parbat we propose a model in which this intense metamorphic and structural reworking of crustal lithosphere is a consequence of strain focusing caused by significant erosion within deep gorges cut by the Indus and Tsangpo as these rivers turn sharply toward the foreland and exit their host syntaxes. The localization of this phenomenon at the terminations of the Himalayan arc owes its origin to both regional and local feedbacks between erosion and tectonics.

Geophysical evidence from the MELT area for compositional controls on oceanic plates

Magnetotelluric and seismic data, collected during the MELT experiment at the southern East Pacific Rise, constrain the distribution of melt beneath this mid-ocean-ridge spreading centre and also the evolution of the oceanic lithosphere during its early cooling history. Here we focus on structures imaged at distances ∼100 to 350 km east of the ridge crest, corresponding to seafloor ages of ∼1.3 to 4.5 million years (Myr), where the seismic and electrical conductivity structure is nearly constant and independent of age. Beginning at a depth of about 60 km, we image a large increase in electrical conductivity and a change from isotropic to transversely anisotropic electrical structure, with higher conductivity in the direction of fast propagation for seismic waves. Conductive cooling models predict structure that increases in depth with age, extending to about 30 km at 4.5 Myr ago. We infer, however, that the structure of young oceanic plates is instead controlled by a decrease in water content above a depth of 60 km induced by the melting process beneath the spreading centre.

3-D resistivity forward modeling and inversion using conjugate gradients

We have developed rapid 3-D dc resistivity forward modeling and inversion algorithms that use conjugate gradient relaxation techniques. In the forward network modeling calculation, an incomplete Cholesky decomposition for preconditioning and sparse matrix routines combine to produce a fast and efficient algorithm (approximately 2 minutes CPU time on a Sun SPARC‐station 2 for 50 × 50 × 20 blocks). The side and bottom boundary conditions are scaled impedance conditions that take into account the local current flow at the boundaries as a result of any configuration of current sources. For the inversion, conjugate gradient relaxation is used to solve the maximum likelihood inverse equations. Since conjugate gradient techniques only require the results of the sensitivity matrix A˜ or its transpose A˜T multiplying a vector, we are able to bypass the actual computation of the sensitivity matrix and the inversion of A˜TA˜, thus greatly decreasing the time needed to do 3-D inversions. We demonstrate 3-D resistivit...

Magnetotelluric evidence for crustal suture zones bounding the Southern Great Valley, California

A geoelectric section inferred from a regional magnetotelluric study across the Coast Ranges, the Great Valley, and the Sierra Nevada reveals significant variations in electrical resistivity. Zones of lower resistivity interpreted at depths from 10 km to at least 30 km lie near mapped geologic boundaries between the Coast Ranges and the Great Valley and beneath the eastern side of the Great Valley. The former boundary is inferred by others to separate the subduction complex of the Coast Ranges from the mafic basement of the Great Valley. The lower resistivities are most likely associated with metasediments trapped between the Coast Ranges ophiolite and the former oceanic crust beneath the Great Valley. The latter boundary is problematic, but may be evidence for a deep metasedimentary section trapped between the ophiolites beneath the Great Valley and the granitic rocks of the Sierra Nevada. The lack of change in the magnetotelluric phase across the Great Valley indicates that a suture zone marked by lower resistivities is unlikely to be present beneath the valley. However, this does not preclude the existence of a resistive suture zone.

Variations in the electrical conductivity of the upper mantle beneath North America and the Pacific Ocean

Variations in the electrical conductivity of the mantle beneath Carty Lake in the Canadian Shield, Tucson in the southwest United States, and Honolulu and Midway in the north central Pacific were determined through the inversion of long-period magnetotelluric and geomagnetic depth sounding data. Inversion of computed response functions is carried out using a minimum structure, regularized approach. The upper mantle beneath Carty Lake is approximately an order of magnitude more resistive than the upper mantle beneath Tucson and nearly 1.5 orders of magnitude more resistive than Honolulu and Midway Island. Inversions were also carried out where the minimum structure constraint was removed at known upper mantle discontinuities. These models show a jump in conductivity of ∼1.5 orders of magnitude across the 660 km discontinuity, a result that is consistent with laboratory experiments on realistic mantle assemblages. Mantle conductivity profiles at Carty Lake are significantly more resistive than those at Tucson, Honolulu, and Midway to depths of ∼300–400 km. These observations likely reflect differing thermal states, the presence (or absence) of partial melt and volatiles, and may also be related to chemical differences between depleted and undepleted upper mantle. The observed conductivity variations may be interpreted as lateral variations in temperature, partial melt, and/or dissolved hydrogen in olivine.

Conjugate direction relaxation solutions for 3-D magnetotelluric modeling

In recent years, there has been a tremendous amount of progress made in three‐dimensional (3-D) magnetotelluric modeling algorithms. Much of this work has been devoted to the integral equation technique (e.g., Hohmann, 1975; Weidelt, 1975; Wannamaker et al., 1984; Wannamaker, 1991). This method has contributed significantly to our understanding of electromagnetic field behavior in 3-D models. However, some of the very earliest work in 3-D modeling concentrated on differential methods (e.g., Jones and Pascoe, 1972; Reddy et al., 1977). It is generally recognized that differential methods are better suited than integral equation methods to model arbitrarily complex geometries, and consequently this area has recently been receiving a great deal of attention (e.g., Madden and Mackie, 1989; Xinghua et al., 1991; Mackie et al., 1993; Smith, 1992, personnal communication). Differential methods lead to large sparse systems of equations to be solved for the unknown field values. It is possible to use relaxation al...

3-D magnetotelluric inversion for resource exploration

Summary The standard approach to solving nonlinear geophysical inverse problems is by iterative, linearized inversion, which when run to completion minimizes an objective function over the space of models. This method, however, requires computing a Jacobian (sensitivity) matrix and solving a nonsparse, linear system on the model space at each inversion iteration. These computational tasks, while tractable for one-dimensional (1-D) and two-dimensional (2-D) inverse problems, are prohibitive for larger and more complicated three-dimensional (3-D) problems. For 3-D magnetotelluric (MT) inversion, we use the method of nonlinear conjugate gradients (NLCG) applied directly to the minimization of the objective function. Given the structure of the MT problem, the NLCG method replaces computation of the Jacobian matrixand solution of a large linear system with computations equivalent to only three forward problems per inversion iteration, dramatically increasing the speed to convergence. The algorithm has been tested on both synthetic and real data. In both cases, the results are in good agreement with either the actual model or with known geology. The computation times are modest, taking approximately 10–12 hours on a 400 MHz desktop computer to invert 100 stations at 5 frequencies using 20 NLCG iterations.