Sharp boundary inversion of 2D magnetotelluric data

Geophysical Prospecting, 1999, 47, 469–486 Sharp boundary inversion of 2D magnetotelluric data1 Torquil Smith,2 Michael Hoversten,2 Erika Gasperikova2 and Frank Morrison2 Abstract We consider 2D earth models consisting of laterally variable layers. Boundaries between layers are described by their depths at a set of nodes and interpolated laterally between nodes. Conductivity within each layer is described by values at a set of nodes fixed within each layer, and is interpolated laterally within each layer. Within the set of possible models of this sort, we iteratively invert magnetotelluric data for models minimizing the lateral roughness of the layer boundaries, and the lateral roughness of conductivities within layers, for a given level of data misfit. This stabilizes the inverse problem and avoids superfluous detail. This approach allows the determination of boundary positions between geological units with sharp discontinuities in properties across boundaries, while sharing the stability features of recent smooth conductivity distribution inversions. We compare sharp boundary inversion results with smooth conductivity distribution inversion results on a numerical example, and on inversion of field data from the Columbia River flood basalts of Washington State. In the synthetic example, where true positions and resistivities are known, sharp boundary inversion results determine both layer boundary locations and layer resistivities accurately. In inversion of Columbia flood basalt data, sharp boundary inversion recovers a model with substantially less internal variation within units, and less ambiguity in both the depth to base of the basalts and depth to resistive basement. Introduction In many instances of interpretation of field data, an investigator may suspect that the geology of a field area lends itself to approximation by a model consisting of a few layers of laterally varying thickness and perhaps laterally varying conductivities. This kind of model is particularly appropriate for many areas of current interest in petroleum exploration world-wide. Magnetotellurics (MT) are most often considered in petroleum exploration in cases which are difficult for seismic imaging. A near-surface unit of high resistivity and velocity (salt, basalt or carbonate) overlying prospective sediments above 1 2 Received February 1998, revision accepted January 1999. Department of Materials Science and Mineral Engineering, University of California, 577 Evans Hall, Berkeley, CA 94720, USA. q 1999 European Association of Geoscientists & Engineers 469 470 T. Smith et al. basement is the most common case. In such cases, the geometry of the base of the resistive unit and of the basement surface are primary interpretational goals. Recent inversions have been successful in finding smoothly varying two-dimensional models fitted to magnetotelluric data (e.g. deGroot-Hedlin and Constable 1990; Smith and Booker 1991). These inversions have been explicitly formulated to minimize some measure of the roughness of a conductivity model for some level of squared data misfit. For example, minimizing  2  2 ∂ ∂ logðjÞ þ logðjÞ dx dz þ bjrj2 ; ∂x ∂z where j(x, z) is the conductivity, jrj2 is the squared data misfit and b > 0 is a trade-off parameter, results in a model which is smooth in both vertical and horizontal directions. For non-linear problems, such as inversion of 2D magnetotelluric data, inversions of this sort are generally made iteratively, linearizing about some current model at each iteration, and a sequence of models is generated. To avoid unnecessary difficulties with non-linearities, inversions are generally started from relatively smooth models and a relatively small value of b. As an inversion progresses, greater emphasis is placed on fitting data by increasing the value of b, resulting in models with progressively more detailed structures. The underlying assumption is that by varying b slowly enough, changes to the conductivity model will be small enough for equations based on linearization to result in progressively ‘better’ models, leading to a model, smoothest for its level of misfit, within a moderate number of iterations. Smooth models fitted to data generated from models with sharp boundaries may not resolve the positions of boundaries as well as desired. A smoothest model adequately fitting a data set makes evident resolution limits of the data set in the sense that the data by themselves certainly do not require sharper transitions than are found in such a model. However, if it is suspected that a field area is comprised of a few relatively uniform geological units with abrupt changes in conductivity between units, then a model parametrized in terms of sharp boundaries between units of differing conductivity may be more appropriate. A number of 2D inversions have been formulated along these lines (e.g. Eysteinsson 1986; Marcuello-Pascual, Kaikkone and Pous 1992). A fairly standard tactic has been to parametrize the location of boundaries, to compute partial derivatives of data with respect to parameters, and to invert iteratively for changes to the parameters using a damped matrix inversion. In such a method, data may not uniquely resolve all desired parameters of an inversion: parameter values may depend on implementation in an arbitrary manner and the resulting models may contain unresolvable details. We minimize unresolvable details by explicitly seeking models that minimize the roughness of the boundaries between units, as well as the roughness of resistivity variations within units. Method In principle, arbitrarily complicated 2D models can be developed by specifying the conductivities of different geological units and specifying the positions of their q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 Sharp boundary inversion 471 boundaries at a series of points along the boundaries. We restrict our attention to models comprised of layers with laterally variable thicknesses, with layer interface positions specified at nn nodes along each interface and interpolated between nodes. For a model consisting of nb layers over a basal half-space, each layer interface is specified by its depth at nn nodes. The series of depths for nb interfaces yields an nb nn long vector m of parameters, ðz11 ; z12 ; . . . ; z1nn ; z21 ; . . . ; znb nn Þt ; ð2aÞ where zij is the depth of the j th node of the i th interface, and superscript t denotes transpose. If conductivities are unknown, nb þ 1 unknown layers and half-space conductivities jj may be appended to m. For layers with laterally varying conductivities, we specify the conductivity within each layer at the nn nodal positions and interpolate it horizontally, and append these additional parameters to m, obtaining
t ÿ ð2bÞ z11 ; . . . ; z1nn ; . . . ; znb nn ; logðj11 Þ; . . . ; logðj1nn Þ; . . . ; logðjnb nn Þ ; where jij is the conductivity of the i th layer at the j th node. Once a model parametrization has been chosen, an inversion is specified by choosing a functional to minimize. For variable thickness layered models with specified layer conductivities, we might choose to minimize a functional such as  2  2 d d z1 ðxÞ dx þ . . . þ znb ðxÞ dx þ bjrj2 ; ð3Þ dx dx for layer boundaries written as continuous functions zi(x) of horizontal position. For models, specified by a finite number of parameters, such as vector (2a), we may approximate derivatives by finite differences and integrals by sums. Defining an nb ? (nn ¹ 1) by nb nn roughening matrix R as 2 6 6 6 6 6 6 6 6 6 6 6 4 ¹1=Dx11 1=Dx11 0 0 ¹1=Dx12 1=Dx12 ... 0 0 ... 0 ... 0 0 0 ... 0 ¹1=Dx1nn ¹1 ... ... 1=Dx1nn ¹1 0 0 ¹1=Dx21 ... ... 0 1=Dx21 ... 0 0 ¹1=Dxnb nn ¹1 1=Dxnb nn ¹1 3 7 7 7 7 7 7 7 7 7 7 7 5 with Dxij being the horizontal spacing between the j th and ( j þ 1)th nodes of the i th interface, then the first nn ¹ 1 rows of the product Rm approximate the derivative dz1/dx along the nn ¹ 1 segments between the nodes of the first interface, the next nn ¹ 1 rows approximate the derivative dz2/dx at the nn ¹ 1 segments along the second interface, and so forth. Weighting the rows of the roughening matrix by the would weight the squares of the elements of Rm by the corresponding Dx1/2 ij corresponding interval lengths Dxij, so that jRmj2 ; mt Rt Rm would approximate q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 472 T. Smith et al. the sum of integrals in derivatives, 2 ¹w11 w11 6 0 ¹w12 6 6 6 6 6 ... Rb ¼ 6 0 6 6 0 6 6 4 0 (3). We use a slightly different weighting of approximate 0 ... w12 0 ... 0 ¹w1nn ¹1 ... 0 ... w1nn ¹1 0 0 ¹w21 ... 0 ... 0 w21 ... ... 0 0 ¹wnb nn ¹1 wnb nn ¹1 3 7 7 7 7 7 7 7; 7 7 7 7 5 ð5Þ and choose weights wij constant for each interface: wi ; wi1 ¼ wi2 ¼ . . . ¼ win n ¹ 1. When laterally varying layer conductivities are included as unknown model parameters, as in parametrization (2b), the roughening matrix is extended with another nb nn columns and nb ? (nn ¹ 1) rows, so that " # Rb 0 R¼ ; ð6Þ 0 Rc where Rc is of the same form as Rb with a possibly different set of weights, w0i, i ¼ 1, . . ., nb. The submatrix Rc multiplies layer conductivity parameters in the product Rm. The lower nb ? (nn ¹ 1) rows of Rm approximate the horizontal derivative d log (j)/dx within the layers. For models with laterally varying layer conductivities, when minimizing jRmj2 þ bjrj2 ; ð7Þ w0i, used in forming Rb and Rc, affect the the relative magnitudes of weights wi and relative importance placed on the smoothness of layer boundaries and layer conductivities, respectively. The relative magnitudes of weights wi, w0i and parameter b affect the relative importance placed on model smoothness and data fit, respectively. Magnetotelluric data depend on model parameters non-linearly, so we use an iterative approach to minimizing (7). We write measurements of apparent resistivity and phase at ns sites and nf frequencies each as an element of an nd long data vector d(obs), and the response of a model m as d. Letting Fi be the matrix of data partial derivatives with respect to model parameters evaluated at some model mi with response di, then linearizing about mi gives d ¹ di ¼ Fi ? ðm ¹ mi Þ; ð8Þ where second-order terms in (m ¹ mi) have been neglected. Writing residuals of the measured data with respect to model response d as a vector, r ; dðobsÞ ¹ d; ð9Þ q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 Sharp boundary inversion 473 and using linearization (8) gives r ¼ d0 ¹ Fi m; ð10Þ where d0 ; dðobsÞ þ Fi mi ¹ di ð11Þ on the i th iteration. We use this expression for r as a function of m in expression (7). The minimum of (7) is a stationary point with respect to perturbations dm: at the minimum the difference between (7) written for m þ dm and r þ dr, where r þ dr ¼ d0 ¹ Fi ? ðm þ dmÞ; ð12Þ and written for m and r, must vanish to the first order in dm, for all possible dm. That is, the coefficient of dm in (7) must vanish. This requirement yields ðb Fit Fi þ Rt RÞm ¼ b Fit d0 : ð13Þ For np model parameters, this is an np by np system of equations. We use this equation to estimate a new model m at each iteration, based on partial derivatives updated about the previous model mi. Additional details are given in the Appendix. For very small values of b, the matrix on the left side of equation (13) approaches Rt R, which is singular for the two examples given above, R ¼ Rb for parametrization (2a) and equation (6) for parametrization (2b). The rows of Rb form a set of nb ? (nn ¹ 1) linearly independent vectors, so Rt R has rank nb ? (nn ¹ 1) in the first example, and twice this in the second. In these examples Rt R is an nb nn or a 2nb nn square matrix, so it has nb or 2nb zero eigenvalues, respectively. The jRmj2 term in the object function (7) penalizes components of m in the directions of eigenvectors of large eigenvalues of Rt R, that is, in the directions of right singular vectors of large singular values of R. Components of m in the directions of eigenvectors of zero eigenvalues of Rt R are undamped. If some of these undamped directions are poorly constrained by the data, as are directions corresponding to eigenvectors of small eigenvalues of Ft F, inversion results may not be particularly satisfying. We augment the roughening matrix R so that Rt R has fewer zero eigenvalues and (Ft F þ Rt R/b) has fewer undamped directions. For interface depth inversion (parametrization (2a)), the roughening matrix R ¼ Rb is insensitive to the means of the depths to each interface, S zij /nn summed over j ¼ 1, nn for i ¼ 1, nb . This suggests augmenting R to add some penalty based on mean interface depths. We augment R with nb ¹ 1 additional rows of the form ð0; . . . ; 0; ¹wdz ; . . . ; ¹wdz ; wdz ; . . . ; wdz ; 0; . . . ; 0Þ; ð14Þ where, in the i th additional row, the negative entries multiply the nn depths to the nodes of the i th interface and the positive entries multiply the nn depths to the nodes of the (i þ 1)th interface in the product Rm, yielding a weighted first difference between average interface depths in the two layers. With this augmentation, Rt R only has one remaining zero eigenvalue, and its eigenvector corresponds to the mean of all interface q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 474 T. Smith et al. depths. This choice of augmentation to Rb introduces a preference for models with the average interface depth in each layer close to those in adjacent layers, that is, a preference for thin layers. Instead of parametrizing the model m in terms of depths to interfaces (2a), we may parametrize it in terms of depth below the previous interface, that is, in terms of the layer thicknesses ðDz11 ; Dz12 ; . . . ; Dznb nn Þt ; ð15Þ where Dzij ; zij ¹ zi ¹ 1 j is the thickness of the i th layer at the j th horizontal node. In this case, the augmentation to Rb introduces a preference for uniform average layer thicknesses. Parametrizing the model m in terms of logarithms of layer thicknesses
t ÿ ð16Þ logðDz11 Þ; logðDz12 Þ; . . . ; logðDznb nn Þ ; allows layer thicknesses to remain positive. Synthetic data example We generated 23 sites of TM-mode data at 26 frequencies logarithmically spaced from 0.001 Hz to 100 Hz, from a simple layered 2D model, with layer boundaries as shown by white lines in Fig. 1. The model consists of a surface layer of 300 Qm (e.g. basalt) over a 20 Qm layer (e.g. sediments) over a 100 Qm basement at 12 km depth. The surface layer thins from 2.5 km to 0.7 km over a lateral distance of 12 km. The model response was calculated using the same forward modelling code (Wannamaker and Stodt 1987) as used in forward modelling steps within our inversion. To simulate measurement errors, 5% Gaussian noise was added to apparent resistivities and impedance phases. In inversion, the model was parametrized in terms of interface depths, log(zij) ( j ¼ 1, 23), (i ¼ 1, 2), at nodes directly below the data sites, and layer and basement conductivities, log(jij), at corresponding nodes ( j ¼ 1, 23), (i ¼ 1, 3). Default roughening matrix weights wi ¼ i and w0i ¼ i were used, placing increased emphasis on smoothness with increasing depth. Default vertical roughening weights of wdz ¼ i/nn /10 were used. Equal emphasis was placed on smoothness of conductivity within layers, and smoothness of layer interfaces (wi ¼ w0i). A goal of one standard error rms misfit was chosen for inversion corresponding to the noise level in the input data, and this was achieved in the final iterations of the inversion. Sharp boundary inversion (SBI) results are shown in Fig. 1a. Layer boundaries are shown in black. The resistivities plotted are those resulting from projection of our variable thickness/variable resistivity layer models on to a fixed finite-element mesh used in calculating their responses (see Appendix). Basalt resistivities (surface layer) are matched within 4% over the length of the profile. The base of this layer is matched to within 150 m throughout the model. Sediment resistivities (second layer) are matched within 5% throughout the entire unit, and basement resistivities are matched within 2%. The particular result shown was started from a 1D model with uniform layer resistivities of 200, 160 and 300 Qm. Starting the SBI inversion from other initial q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 Sharp boundary inversion 475 Figure 1. (a) Sharp boundary inversion of data from a three-‘layer’ model with sloping interface. Black lines: inversion layer interfaces. White lines: true layer interfaces. (b) Minimum structure smooth inversion of same data. White lines: position of layer interfaces in true model. models has given almost identical results, giving some confidence that the algorithm has found a global minimum of the object function. However, in inversions such as this, in which model geometry is included in the unknown parameters, we expect that local minima may be a greater concern than they seem to be in inversions in which model geometry is held fixed. For comparison, the same synthetic data have been inverted using a smoothly varying inversion algorithm (deGroot-Hedlin and Constable 1990), started from a 300 Qm half-space and resulting in the model shown in Fig. 1b. In addition to the deGroot-Hedlin-Constable (1990) algorithm, we have made extensive use of the RRI algorithm of Smith and Booker (1991) and a non-linear conjugate gradient algorithm of Mackie et al. (1997), both of which also invert for smooth conductivity distributions. In our experience, there are three main observations that can be made about smooth inversions of data from models with sharp boundaries, which are illustrated to varying q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 476 T. Smith et al. degrees in Fig. 1. Firstly, in smooth inversions of data from models with interfaces that vary considerably in depth, different contour intervals correspond to the actual interface position at different depths: the deeper the interface, the lower the resistivity of the corresponding contour. For example, on the right in Fig. 1b, the upper interface lies at 720 m depth, with resistivity between 250 and 350 Qm, while on the left it lies at 2550 m depth with resistivity between 30 and 40 Qm, in the smooth inversion result. Secondly, smooth inversion of data from models with sharp boundaries tends to overshoot resistivities on either side of a boundary. In Fig. 1b, the smooth inversion results overshoot by about 20%, attaining highs of about 350 Qm for the 300 Qm basalt, and lows of about 16 Qm for the 20 Qm sediments. Thirdly, smooth inversion results generally have increased smoothing with increasing depth. In Fig. 1b, the transition from sediment to basement resistivities is more spread out than the transition from basalt to sediment resistivities in the smooth inversion results. The data by themselves are unable to resolve whether the transition to basement resistivities is smooth or abrupt. However, assuming that the transition is abrupt, as we do in Fig. 1a, the depth and resistivity of basement are quite accurately recovered. Field data example We consider magnetotelluric data from a survey line somewhat east of Ellensburg, Washington, shown in juxtaposition with topography and interpreted structure in Fig. 2 (line 400). The line consisted of 34 sites, spaced 242 m apart with in-line electric field dipoles laid end to end. Original processing frequencies differ from site to site, so the data were binned to 5 frequencies per decade yielding 25 frequencies common to all sites. Line 400 was laid out as a dip line to be approximately perpendicular to strike. Geoelectric strike interpreted from polar diagrams (not shown) varies somewhat along the line, but is generally perpendicular to the line direction. We interpret the Zxy impedances, with E in-line and B perpendicular, as TM-mode data, and it is plotted in terms of apparent resistivity and phase in Fig. 3. Interpreting this polarization as 2D TM-mode data, the obvious site-to-site shifts in apparent resistivity magnitude can be fitted by lateral resistivity variations at the earth’s surface in a 2D model. Since the data are endto-end dipole data, we do not expect them to be affected by surface structures much smaller than the site spacing. For want of rigorous impedance variance estimates, we have assumed 10% uncertainty in both apparent resistivity and phase data. In the vicinity of the survey, approximately two kilometres of resistive flood basalts overlie several kilometres of more conductive sediments, in turn overlying a more resistive basement. More details of the geological setting and the MT data set are given by Morrison et al. (1996). The depth to the sediments and their thickness are of particular concern as gas deposits have been found and exploited within the sediments at other locations in the Columbia Basin. A two-layer-over-basement model (as in Fig. 1) includes the parameters of greatest interest. We consider the site-to-site shifts observed in the line 400 magnitude data to be most probably due to lateral variations at the earth’s surface, so add a near-surface layer to allow this to be modelled. We q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 Sharp boundary inversion 477 Figure 2. Location of Columbia basalt data line 400, with respect to topography (greyscale) and anticlinal structures interpreted from surface geology (solid lines), from Morrison et al. (1996). also add a deep conductor to enable us to fit the low-frequency roll-off of apparent resistivities. The SBI algorithm was started from a four-layer-plus-basement model with interface depths of 0.2, 2, 4 and 11 km. Starting layer resistivities were 300, 300, 20, 150 and 20 Qm. Parameter nodes were below the measurement sites. This yields a total of 136 interface nodes and 170 resistivity nodes. The final SBI inverse model which reached an rms misfit of 2.0 standard errors is shown in Fig. 4. The upper two layers comprise the basalt flow sequence with the upper part more conductive (< 100 Qm) than the lower part (< 300 Qm). The increased conductivity is consistent with the expected effects of weathering near the surface. The sediments are represented by the third layer with resistivities ranging from 3 to 28 Qm. Basement resistivities range from 80 to 300 Qm. The deep conductor (not shown) averaged 5 Qm with an average depth of 15 km. Inversion model response and input data are plotted together in x ¹ y plots site by site in Figs 5 and 6. For comparison, we inverted the same data set using the deGroot-Hedlin and Constable (1990) smooth conductivity distribution inversion algorithm. Starting from q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 478 T. Smith et al. q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 Figure 3. Columbia basalt line 400, Zxy impedance data, with Ex in-line direction. (a) Apparent resistivity. (b) Phase. Sharp boundary inversion 479 Figure 4. Sharp boundary inversion of data from Columbia basalt line 400. White lines: inversion interface location. Colour: resistivity of sharp boundary layered model projected on to finite-element mesh. a 300 Qm half-space, the algorithm resulted in a model with an rms misfit of 1.9 standard errors, shown in Fig. 7. We have superimposed the SBI inverse layer boundary positions. Differences between the two models underline the difficulties involved in making a discrete geological unit interpretation of a smooth conductivity distribution inversion. The deGroot-Hedlin-Constable (1990) smooth inversion shows considerably more lateral variation in resistivities, both within the basalts and within the underlying sediments. One might interpret a somewhat deeper basalt base from the smooth inversion than is evident in the sharp boundary inversion results (bottom of second layer). Basement relief is quite consistent between the two inversions, with the smooth inversion yielding lower resistivities at a given depth than the sharp boundary inversion, consistent with the synthetic results in the previous section. Conclusion If relatively little is known about a survey area, interpretation in terms of a minimum structure smooth inversion such as those of deGroot-Hedlin and Constable (1990) or Smith and Booker (1991) quickly indicates what can be inferred about the study area from a data set by itself. The blurriness of the results of a minimum structure smooth inversion gives some idea as to resolution limits of a given data set. When it is suspected that the field area is made up of a few relatively homogeneous units, the positions of boundaries can be recovered more precisely using a sharp boundary inversion. Inverting directly for boundary positions in a sharp boundary inversion makes interpretation for structural geology easier. q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 480 T. Smith et al. q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 Figure 5. TM-mode apparent resistivity response of sharp boundary inversion (Fig. 4) as a function of frequency at individual measurement sites, from top to bottom, left to right. Smooth curves: model response. Erratic curves: Columbia basalt line 400 data. Sharp boundary inversion 481 q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 Figure 6. TM-mode impedance phase response of sharp boundary inversion (Fig. 4) as a function of frequency at measurement sites, from top to bottom, left to right. Smooth curves: model response. Erratic curves: Columbia basalt line 400 data. 482 T. Smith et al. Figure 7. Minimum structure smooth inversion of Columbia basalt line 400 data. White lines: sharp boundary inversion interface positions from Fig. 4. A sharp boundary inversion provides MT interpreters with a tool specifically designed for structural interpretation. By parametrizing inverse models in terms of boundaries between units which may possess large contrasts in conductivity, inversions which unambiguously locate such boundaries can be produced. In the case of petroleum exploration, where the base of a sequence of relatively uniform basalt, salt or carbonate is the goal of an MT survey, a sharp boundary inversion provides advantages over a smooth inverse model which leaves the interpreter with the need to interpret the smooth inverse model for the location of structural interfaces. Acknowledgements The Engineering Geosciences Program at UCB is grateful to the Marine Magnetotelluric Consortium of AGIP, BHP, BP, British Gas, Chevron, Texaco and Western Atlas, in cooperation with the US Department of Energy, Offices of Energy Research and of Computational and Technology Research, Laboratory Technology Applications Division under Contract DE-AC03-76SF0098, for their support of research. We thank Shell Exploration for the Columbia River basalt data used in this paper. Appendix A non-linear minimum roughness inversion algorithm We approach the minimum of expression (7) iteratively, starting at some model mi computing the partial derivative matrix Fi numerically at mi, and solving equation (13) for an estimated new model m̂i þ 1 ¼ m. Equation (13) is based on linearization (8), q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 Sharp boundary inversion 483 which assumes that the change in the model (m̂i þ 1 ¹ mi) is small enough to neglect second-order terms in (m̂i þ 1 ¹ mi). We consider an iteration successful if, on computing the response di þ 1 of a model mi þ 1, the objective function (7) decreases for the current value of b, i.e. if jRmiþ1 j2 þ bjriþ1 j2 # jRmi j2 þ bjri j2 : ðA1Þ As the model where partial derivatives are evaluated, mi, may not yet be at a minimum of expression (7) for any value of b, adjusting b and recomputing m̂i þ 1 does not guarantee that (m̂i þ 1 ¹ mi) can be made small enough for the condition (A1) to hold for mi þ 1 ¼ m̂i þ 1. Taking a partial step in the desired direction by letting miþ1 ¼ mi þ a ? ½m̂iþ1 ¹ mi ÿ; ðA2Þ with 0 < a # 1, we can choose a small enough so that condition (A1) can be made to hold for mi þ 1, assuming that equation (13) is non-singular. Smith and Booker (1988) give a simple criterion for choosing between reducing a and reducing b to reduce the step size. In their formulation, the squared misfit can be written as a function of the trade-off parameter b with the form 2 nd  X gk 2 ; ðA3Þ jr̂iþ1 ðbÞj ¼ 1 þ blk k¼1 (their equation A9), where lk $ 0 and gk are known, for the estimated model m̂i þ 1(b) considered as a function of b, based on a linearization similar to equation (8). They solve an equation of the form (A3) for the value bs which linearization predicts should result in the smoothest model m̂i þ 1(bs) with the current level of misfit jr̂i þ 1(bs)j2 ¼ jrij2. They divide the model change a ? [m̂i þ 1 ¹ mi] into two parts, a ? [m̂i þ 1(b) ¹ m̂i þ 1(bs)] and a ? [m̂i þ 1(bs) ¹ mi]. The first part is proportional to the difference between an estimated model m̂i þ 1(b), with b chosen so as to reduce the misfit, and the model estimated to be the smoothest model with the current misfit, so representing a model improvement step. The second part represents a model smoothing step as it attempts to smooth the current model without reducing its misfit. When a step size needs to be reduced to lower the object function then, if jm̂i þ 1(b) ¹ m̂i þ 1(bs)j > jm̂i þ 1(bs) ¹ mij, the improvement part of the step is diminished by choosing b closer to bs, otherwise the entire step size is diminished by reducing a. We use the same criterion here. Smith and Booker (1988) derived their equation of form (A3) for a 1D continuous profile inverse problem. Smith (1988) derived the analogous equation for a minimum roughness discrete inverse such as is presented here. Calculation of coefficients gk and lk involves finding the eigenvectors and eigenvalues of an nd by nd matrix computed from Fi and R. An alternative method to solve for b giving a chosen level of the estimated squared misfit jr̂i þ 1(b)j2 is to solve equation (13) for model estimates m̂i þ 1(b) for a few values of b using Cholesky decompositions of (bFti Fi þ Rt R), use equation (10) to evaluate their estimated residuals r̂i þ 1(b), and interpolate between q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 484 T. Smith et al. squared misfit values jr̂i þ 1(b)j2 using an interpolation function based on a single term of the form of expansion (A3):  2 g 2 ; ðA4Þ jr̂ðbÞj < 1 þ bl with g and l chosen to match the calculated jr̂i þ 1(b)j2 at two values of b. Equation (A4) can be solved for b for desired levels of estimated misfit. Letting b i and b is be two values of b used in the previous iteration, m̂i þ 1(b i) and m̂i þ 1(b is) be the corresponding models estimated using equation (13), and r̂i þ 1(b i) and r̂i þ 1(b is) be their respective estimated residuals from equation (10), we interpolate jr̂i þ 1(b)j2 using (A4), and solve for b si þ 1 yielding an estimated squared misfit jr̂i þ 1(b is þ 1)j2 equal to the current level of squared misfit jri j2. We also use the interpolating function (A4) to solve for the reciprocal damping factor b gi þ 1 which yields an estimated residual jr̂i þ 1(b igþ 1)j2 equal to the level eventually desired. The reciprocal damping factor b i þ 1 used in an iteration is chosen between these values. In particular, we choose b i þ 1 based on the scales of Rt R and Fti Fi : biþ1 ¼ z ? traceðRt RÞ=traceðFti Fi Þ; ðA5Þ where the trace of a matrix is the sum of its diagonal elements, and z is an adjustable parameter. In general, we start with z ¼ 1 on the first iteration, and z is later adjusted in tandem with a to control the step size, eventually attaining values giving b i þ 1 ¼ b igþ 1 on later iterations. In forming equations (11) and (13) at each iteration, we need to compute the partial derivative matrix Fi. We do this by finite differencing. We use subroutines based on the finite-element method of Wannamaker and Stodt (1987) to solve the electromagnetic induction forward problem on a fixed mesh. Our variable thickness layered models m are projected on to the fixed mesh of the finite-element grid, and element conductivities are assigned in proportion to element portions within each layer. These subroutines set up linear systems of equations of the form Ax ¼ b ðA6Þ for each frequency and mode modelled, and we factorize the coefficient matrix A as A ¼ LU where L is lower triangular and U is upper triangular, to solve Ax ¼ b, by solving Lx0 ¼ b and Ux ¼ x0 . In computing partial derivatives we use the distributive law to write perturbed systems as A ? ðx þ dxÞ ¼ ðb þ dbÞ ¹ dA ? ðx þ dxÞ; ðA7Þ A ? ðx þ dxÞ ¼ ðb þ dbÞ þ b ¹ ðA þ dAÞx ¹ dA dx; ðA8Þ or where (A þ dA) and (b þ db) are calculated for a perturbed model (mi þ dm), vectors b and x are obtained from previous calculation of the unperturbed case, and we solve (A8) for (x þ dx), neglecting the second-order term dA dx. This allows (A8) to be q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 Sharp boundary inversion 485 solved using the LU factors of A previously calculated in solving the unperturbed forward problem. Even with reuse of LU factors, computing partial derivatives takes a significant amount of time. To save computation, we approximate some partial derivatives with their values determined in previous iterations (evaluated at previous models). To avoid the effects of a systematic ordering of parameters in the model vector m, we evaluate the partial derivatives cycling through a set of np random orthogonal directions. Letting ek, (k ¼ 1, np) be unit vectors in np orthogonal directions, and Jek be the vector of partial derivatives of the data d computed for a small model perturbation dm in direction ek, then F 0 ¼ F þ ðJek ¹ Fek Þ etk ðA9Þ updates the partial derivative matrix in this direction. Cycling through √ all np directions — results in updating the entire matrix. We update F in approximately np orthogonal directions per iteration. In addition to these, we also update F in the directions of the most recent improvement step (mi(bi) ¹ mi(b is)), and the orthogonal component of the corresponding smoothing step (mi(b is) ¹ mi ¹ 1). In short, after updating the partial derivative matrix, we solve equation (13) for the estimated models m̂i þ 1(b is) and m̂i þ 1(b i) for use in solving for b si þ 1 and b igþ 1, and then solve for m̂i þ 1(b is þ 1) and m̂i þ 1(b i þ 1). Forming the new model mi þ 1 through equation (A2), we solve the full 2D forward problem for the new model mi þ 1 to compute its response di þ 1 and residual ri þ 1. If the new model mi þ 1 reduces the object function (i.e. satisfies (A1)), then it is accepted as the starting point for the next iteration, otherwise the step size is reduced through a or z and F is updated in more directions before computing mi þ 1 and di þ 1 again. If two successive iterations are successful without decreases in a or z, then a or z is increased to increase the size of the next attempted step. References Eysteinsson H. 1986. Two-dimensional inversion of magnetotelluric data. 8th workshop on electromagnetic induction in the Earth and Moon, IAGA working group I-3, Neuchatel, Switzerland. Expanded Abstracts. deGroot-Hedlin C. and Constable S.C. 1990. Occam’s inversion to generate smooth, twodimensional models from magnetotelluric data. Geophysics 55, 1613–1624. Mackie R.L., Livelybrooks D.W., Maden T. and Larsen J.C. 1997. A magnetotelluric investigation of the San Andreas fault at Carrizo Plain, California. Geophysical Research Letters 24, 1847–1850. Marcuello-Pascual A., Kaikkone P. and Pous J. 1992. Two-dimensional inversion of magnetotelluric data with a variable model geometry. Geophysical Journal International 110, 297–304. Morrison H.F., Shoham Y., Hoversten G.M. and Torres-Verdin C. 1996. Electromagnetic mapping of electrical conductivity beneath the Columbia basalts. Geophysical Prospecting 44, 963–986. q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486 486 T. Smith et al. Smith J.T. 1988. Rapid inversion of multi-dimensional magnetotelluric data. PhD dissertation, University of Washington, Seattle. Smith J.T. and Booker J.R. 1988. Magnetotelluric inversion for minimum structure. Geophysics 53, 1565–1576. Smith J.T. and Booker J.R. 1991. Rapid inversion of two- and three-dimensional magnetotelluric data. Journal of Geophysical Research 96, 3905–3922. Wannamaker P.E. and Stodt J.A. 1987. A stable finite element solution for two-dimensional magnetotelluric modelling. Geophysical Journal of the Royal Astronomical Society 88, 277–296. q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 469–486