Magnetotelluric inversion for anisotropic conductivities in layered media

Physics of the Earth and Planetary Interiors 158 (2006) 139–158

Magnetotelluric inversion for anisotropic conductivities in layered media Josef Pek a,∗ , Fernando A.M. Santos b a

Geophysical Institute, Academy of Sciences of the Czech Republic, Boˇcn´ı II/1401, CZ-14131 Prague 4-Spoˇrilov, Czech Republic b Physics Department-CGUL, University of Lisbon, Campo Grande, C8, Piso 6, 1749-016 Lisbon, Portugal Received 3 August 2005; received in revised form 7 December 2005; accepted 15 March 2006

Abstract Electrical anisotropy in the Earth’s crust and upper mantle has recently gained attention as a significant linking factor between electrical models and underlying structural and tectonic patterns. This interest has also motivated new methodological studies into the modelling and inversion for electrically anisotropic structures. We present an algorithm for the inversion of magnetotelluric data over layered anisotropic conductors which is a straightforward extension of the standard Occam 1-D inversion to anisotropic models. Owing to the essential limitation of magnetotellurics to resolve the complete conductivity tensor, we formulate the inversion for azimuthal anisotropy only. We treat the non-linear inverse problem as a multi-criterion minimization of the structure complexity, data misfit and anisotropy. To constrain the structure complexity, we employ the standard roughness penalty as well as non-quadratic penalties of the total variation and gradient support type that produce more focused model sections and thus conform better to the idea about sharp, non-diffuse boundaries of anisotropic structures in the Earth. Application of the anisotropy penalty is crucial for suppressing spurious anisotropy in the inverse models. We use a 2-D extension of the heuristic L-curve method to estimate the quasi-optimal penalty weights. With two non-linear iteration solvers, specifically the reweighted conjugate gradient method and the lagged diffusivity iteration, we can arrive at the minimum of the target functional, for one selected pair of regularization weights, typically after a few tens of iteration steps. To demonstrate the inverse solution, we present two simple yet not completely trivial synthetic examples, the first one based on data generated by a model with two anisotropic layers with discordant strikes, and the other showing possible misinterpretations in case a 1-D inversion with anisotropy is formally applied to data produced by simple 2-D block structures. Field data from the Variscan structures in southern Portugal bear indications, both geological and geophysical, on possible crustal anisotropy due to shearing and graphitization. We present results of a 1-D anisotropic imaging, and a preliminary 2-D modelling, for a short section of five sites over an old suture between the South Portuguese Zone and Ossa Morena Zone, aimed at explaining some essentially non-2-D effects in the data, as, e.g., systematic frequency variations of the regional strike or phases rolling out of their natural quadrants. © 2006 Elsevier B.V. All rights reserved. Keywords: Magnetotelluric method; Anisotropic conductor; Layered medium; Inverse problem; Regularization; Occam inversion

1. Introduction

∗

Corresponding author. Fax: +420 272761549. E-mail addresses: [email protected] (J. Pek), [email protected] (F.A.M. Santos). 0031-9201/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.pepi.2006.03.023

Recent progress in geoelectrical induction studies has demonstrated that electrical anisotropy in the deeper parts of the Earth can be a real and significant factor of the Earth’s internal structure. Though arguments

140

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

are still raging as regards a possibility of distinguishing between effects of anisotropy and manifestations of lateral conductivity variations or local distortions, several well-grounded studies have proved the electrical anisotropy to be fully justified within specific geological and tectonic settings (e.g., Mareschal et al., 1995; Eisel and Haak, 1999). Moreover, those studies have also proven that electrical anisotropy creates a significant link between the geoelectrical models and geodynamic interpretations (for a recent review, see Wannamaker, 2005). Detailed studies of electrically anisotropic Earth structures have induced further interest to new interpretation techniques for anisotropic conductors. As compared with isotropic models, inversion for anisotropic electrical conductivities has to cope with a larger number of variables as well as with more complex interrelations between the parameters. Another important fact is that the anisotropic structures within the Earth are often related to spatially strictly demarcated tectonic zones, and do not conform the general smoothness concept as employed in most of standard Occam inversions based on penalizing model roughness or curvature. First attempts to invert magnetotelluric (MT) data for 1-D anisotropic conductivities were made in 1970s by Abramovici and Shoham (1977) who based their inversion procedure on the generalized inversion technique. More recently, Regis and Rijo (1997, 2000) have used the same inversion procedure, and demonstrated that inverting the electrical data jointly with additional prior information from well logging and other geophysical and geological interpretations is crucial for stabilizing the anisotropic models. Lately, a version of an efficient nonlinear inversion algorithm for a 1-D CSAMT problem with azimuthal anisotropy has been developed and successfully applied to data from a shallow fractured crustal environment by Li et al. (2000). As the non-uniqueness of the inverse problem solution is generally much more serious in the case of anisotropic structures, direct or random search nonderivative optimization procedures may be used with advantage to perform an exhaustive mapping of the whole parameter space. With this aim, Santos and Mendes-Victor (1997) used the simulated annealing algorithm to search for the optimal solution to the inverse problem as well as for characterizing whole classes of data-compatible layered models, with both the azimuthal and dipping anisotropy considered within the layers. Later, the question of ambiguity and equivalency of the conductivity parameters in layered anisotropic media was addressed theoretically by Yin (2003) from the point of view of the inverse problem solution.

In this paper, we concentrate on a non-linear inversion of 1-D MT anisotropic data based on an iterative linearized optimization for electrical parameter estimation in a finely discretized layered Earth. We extend the standard Occam 1-D inverse technique (Constable et al., 1987) to anisotropic layered conductors. In Section 2, we formulate our 1-D MT anisotropic model and summarize principles of the solution to the direct problem and calculation of its parametric sensitivities. Section 3 then deals with the solution to the inverse problem. We study the effect of different stabilizing functionals, both quadratic and non-smooth, for the regularization of the inverse problem and discuss the role of the anisotropy penalization in mitigating excessive anisotropy in inverse models. We further describe the technical procedures of the inversion based on the reweighted conjugate gradient and lagged diffusivity minimization routines. We demonstrate the performance of the inverse algorithm on a synthetic model of a double-strike anisotropic Earth. Another synthetic example is then used to demonstrate possible misinterpretations that can arise from formally imaging simple laterally non-uniform 2-D structures by 1-D anisotropic models. Finally, practical MT curves from MT measurements across Variscan structures of southern Portugal are presented and used to illustrate additional aspects of the anisotropic inversion of real experimental data. 2. Model, forward solution, parametric sensitivities 2.1. Model of the 1-D MT problem for anisotropic conductors Throughout this paper, we will assume a simple Cagniard-Tichonov 1-D MT model extended to anisotropic layered media. The conductive structure (Earth) consists of homogeneous horizontal layers with thicknesses hl , l = 1, . . ., N. The stack of layers is underlain by a homogeneous conductive halfspace. The electrical conductivity in each of the layers, as well as in the basement, is given by a conductivity tensor σ l , l = 1, . . ., N, N + 1. We assume that σ is a symmetric and positive definite tensor anywhere within the Earth. The halfspace above the Earth is filled with perfectly insulating air, i.e., σ 0 = 0. The Earth’s surface defines the xy plane of a cartesian coordinate system, the z axis is directed down into the conductive medium. The electromagnetic field is excited by a primary electromagnetic plane wave that originates from sources at z = −∞ and propagates perpendicularly towards the Earth’s surface.

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

Within each layer, the conductivity tensor σ can be, due to its symmetry, diagonalized and expressed via three principal conductivities, σ 1 , σ 2 , σ 3 , and a rotation matrix, which can be, in turn, decomposed into three elementary Euler’s rotations, successively through angles αL , αD , and αS , ⎛

σxx

σxy

σxz

⎞

⎜ ⎟ σ = ⎝ σyx σyy σyz ⎠ = Rz (−αS )Rx (−αD )Rz (−αL ) σzx σzy σzz ⎞ ⎛ σ1 0 0 ⎟ ⎜ × ⎝ 0 σ2 0 ⎠ Rz (αL )Rx (αD )Rz (αS ), 0 0 σ3 where Rx and Rz are elementary rotation matrices around the coordinate axis x and z, respectively. According to their physical relation to the anisotropy, angles αS , αD , αL can be identified with the anisotropy strike, dip, and slant, respectively. Considering the symmetry of the above MT model, i.e. ∂/∂x = ∂/∂y ≡ 0, Maxwell’s equations in each of the homogeneous subdomains of the model for a frequency ω reduce to ∂Ex ∂z ∂Hy ∂z ∂Ey ∂z ∂Hx ∂z

= iωµ0 Hy , = −Jx = −σxx Ex − σxy Ey − σxz Ez , (2.1) = −iωµ0 Hx , = Jy = σyx Ex + σyy Ey + σyz Ez .

141

where σxz σzx , σzz σyz σzx = σyx − , σzz

σxz σzy , σzz σyz σzy = σyy − , σzz

Axx = σxx −

Axy = σxy −

Ayx

Ayy

(2.4)

with clearly Axy = Ayx for a symmetric conductivity tensor σ. From (2.1) and the system (2.3), we can conclude that the MT field in a layered anisotropic medium depends on the elements of the conductivity tensor through the aggregate conductivities Axx , Ayy , and Axy only. Whatever the particular form of the conductivity tensor σ, the electromagnetic field does not change provided the elements of the 2 × 2 matrix A remain unchanged. Consequently, without any additional information available, the MT field of a plane wave does not allow us to reconstruct the full conductivity tensor in a 1D medium. Only the elements of A can be resolved, which thus represent an equivalent effective conductivity tensor attributed to the individual layers of the model. It can be easily shown that det A = det σ/σ zz > 0. Consequently, A is a symmetric and positive definite 2 × 2 matrix, which can be again factorized in terms of its principal elements, say A1 and A2 , and an elementary rotation, by an effective strike angle βS , around the z coordinate axis,      Axx Axy cos βS − sin βS A1 0 = 0 A2 Ayx Ayy sin βS cos βS   cos βS sin βS × , − sin βS cos βS

The last pair of governing field equations degenerates to which gives Hz = 0,

Jz = σzx Ex + σzy Ey + σzz Ez = 0, (2.2)

expressing simply the absence of a vertical magnetic field and of vertical electric currents anywhere in an anisotropic layered medium. Eliminating, e.g., the magnetic field components from (2.1), we easily arrive at a system of coupled secondorder differential equations for the electric field, ∂ 2 Ex + iωµ0 (Axx Ex + Axy Ey ) = 0, ∂z2 ∂ 2 Ey + iωµ0 (Ayx Ex + Ayy Ey ) = 0, ∂z2

(2.3)

Axx = A1 cos2 βS + A2 sin2 βS , Ayy = A1 sin2 βS + A2 cos2 βS , Axy = Ayx = (A1 − A2 ) sin βS cos βS . Summarizing the above steps, we conclude that the 1-D MT problem for a generally anisotropic layered medium can be always re-formulated as a simpler, but equivalent, problem for an azimuthally (horizontally) anisotropic structure with the horizontal conductivity tensor A defined by (2.4). Any changes in the full conductivity tensor σ that do not affect the elements of A cannot be recognized by MT soundings with a plane wave source field.

142

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

2.2. Principles of the forward solution and evaluation of the parametric sensitivities In the past 30 years, the solution to the 1-D direct MT problem for anisotropic structures has been presented in various forms by a number of authors (e.g., O’Brien and Morrison, 1967; Reddy and Rankin, 1971; Loewenthal and Landisman, 1973; Abramovici, 1974; Dekker and Hastie, 1980). Therefore, we will only briefly summarize the principal points of the direct solution here. For a detailed description of the problem, both theoretical and algorithmic, we refer the reader to Pek and Santos (2002). The system (2.3) has the form of coupled pendula equations and its general solution in any of the homogeneous layers can be expressed as a linear combination of four elementary ‘waves’, exp(±k1,2 z), where the wave 2 = −iωµ A numbers are k1,2 0 1,2 in terms of the principal effective conductivities A1 and A2 . Thus, two pairs of downgoing and upgoing ‘waves’ exist in a generally anisotropic layer, the first one, corresponding to the wave number k1 , ‘slow’, and the other one ‘fast’ (k2 ), provided A1 ≥ A2 . A standard matrix propagation method, based on fitting the fields on the layer boundaries, can now be used to propagate the fields between different depths within the model. For a particular z within an lth layer, the matrix propagation procedure provides a general

fields in a 1-D anisotropic layered medium (e.g., O’Brien and Morrison, 1967; Loewenthal and Landisman, 1973). The matrix expression (2.5) can be further used to derive generalized recurrent formulae for propagating the MT impedances between the boundaries of the anisotropic layers. Considering in particular the top and the bottom of the lth layer, z = zl−1 and z = zl , respectively, we easily have a relation between the respective horizontal fields, Fl−1 (zl , ω) = Sl (hl , ω)Fl (zl , ω). After introducing the impedance tensor by Eh (z, ω) = Z(z, ω)Hh (z, ω), we obtain a system of matrix equations EH Zl−1 Hh,l−1 = (SEE l Zl + Sl )Hh,l , HH Hh,l−1 = (SHE l Zl + Sl )Hh,l ,

with SEE through SHH being 2 × 2 logical sub-blocks l l of the matrix Sl which relate the respective electric and magnetic sub-vectors of the horizontal field between the layer boundaries. Finally, −1

EH HE HH Zl−1 = (SEE l Zl + Sl )(Sl Zl + Sl )

,

(2.6)

where Zl ≡ Z(zl , ω) and Hh,l = Hh (zl , ω). When used recursively, Eq. (2.6) can be used to propagate the impedance tensor from the top of the homogeneous basement through all layers up to the Earth’s surface. The starting tensor on the top of the homogeneous anisotropic halfspace is given by   −(ζ1 − ζ2 ) sin 2βS (ζ1 + ζ2 ) + (ζ1 − ζ2 ) cos 2βS 1 ZN = − , 2 −(ζ1 + ζ2 ) + (ζ1 − ζ2 ) cos 2βS (ζ1 − ζ2 ) sin 2βS

formula for the horizontal field components (e.g., Pek and Santos, 2002), Fl (z, ω) = Sl (zl − z, ω)

N

Sj (hj , ω)DN+1 ,

j=l+1

z ∈ (zl−1 , zl ),

hj = zj − zj−1 ,

(2.5)

where FT = (ETh , HTh ) = (Ex , Ey , Hx , Hy ) is the horizontal field vector, and each of the 4 × 4 matrices Sj (hj , ω) depends on the parameters of the jth layer only. The 4-vector DN+1 contains only two non-zero entries that are arbitrary constants corresponding to the complex amplitudes of the downgoing solution modes in the infinite homogeneous basement. For a particular problem, these constants can be determined, e.g., by normalizing two of the four field components on the Earth’s surface. Eq. (2.5) allows us to propagate the electromagnetic field from the basement of the model through its layers towards the Earth’s surface. It represents the mathematical basis of the matrix propagation procedure for the

with ζ 1 , ζ 2 , βS being the halfspace parameters, ␨m ≡ ζ m,N+1 = km,N+1 /Am,N+1 , m = 1,2, and βS = βS,N+1 . Eq. (2.6) is a basis for the impedance propagation method in a 1-D anisotropic layered conductor (e.g., Abramovici, 1974; Dekker and Hastie, 1980). As, in general, impedances behave more regularly than the fields themselves, the impedance propagation method is preferred in the numerical evaluation of the impedance tensor. Parametric sensitivities of the surface impedances with respect to the parameters of the individual layers are given by partial derivatives ∂Z(z0 , ω)/∂pk , where pk is one of the model parameters of the kth layer. In MT inverse applications for anisotropic 1-D conductivities/resistivities without any specific prior information on the anisotropy, it is reasonable to choose pk ∈ {A1,k , A2,k , βS,k , hk }, which is a quadruplet of independent parameters that completely characterize any anisotropic layer within a stratified medium. The parametric sensitivities of the MT impedance for an anisotropic layered halfspace can be computed directly by differentiating the formulae (2.6) with respect

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

to pk . Since (i) the matrix Sl , and consequently all its subblocks depend solely on the parameters of the current, i.e. the lth layer, and (ii) the MT impedance at the bottom of the lth layer, Zl , does not depend on any parameter of any layer above that boundary, the evaluation of the sensitivities can be very easily integrated into the algorithm for the impedance propagation (Pek and Santos, 2002). As the electromagnetic field vector does not meet the second condition above, sensitivity computations via a field propagation procedure would be more involved. 3. Inversion 3.1. Target function Our approach to the MT inversion for anisotropic layered models is based on the standard Occam idea (e.g., Constable et al., 1987) of minimizing the structure with the data fit used as a constraint in the minimization problem. The parametric functional to be minimized can be written in a general form Φ(m, λ) = Φstruct (m, mprior ) + λ−1 Φdat (m)

2



→ min, with Φdat (m) = WD [dmod (m) − dobs ] , (3.1) where m is the vector of the model parameters, mprior contains prior estimates of the model parameters, dobs are the measured data, WD the covariance matrix of the experimental data, and dmod (m) is the forward solution of the MT problem for the model parameters m. The parameter λ is a regularization parameter which specifies the weight given to the structural constraint Φstruct (m, mprior ) at the expense of the data fit Φdat (m). In our experiments, the parameterization of the model is based on dividing the 1-D section into a logarithmically uniform system of homogeneous layers between two depth limits, which are derived from estimates of the minimum and maximum apparent penetration depths of the electromagnetic field in the medium. Then, the parameters of the model (vector m) are logarithms of the minimum and maximum horizontal resistivities and the direction of the minimum resistivity (anisotropy strike) within each of the layers. We have proven in Section 2.1 that these parameters describe the anisotropic layer completely for MT considerations. For a fine enough sampling of the halfspace, the layer geometry (layer thicknesses) need not be included into the set of the variable parameters. In resolving the parameters of an anisotropic conductor, the secondary (diagonal) impedances play a cru-

143

cial role. Therefore, the experimetal data vector dobs is formed of all the impedance elements rather than the commonly used apparent resistivity and phase data derived from the main (antidiagonal) impedance elements. The aim of introducing the structural penalty Φstruct (m, mprior ) is to prevent the model from showing unnecessary structural features that result from overfitting the experimental data, i.e., from fitting the noise component in the observations. The minimum norm (e.g., Tikhonov and Arsenin, 1977), minimum roughness and minimum curvature (Constable et al., 1987) are the most popular structural penalties used in geoelectrical inversions, both due to their interpretation efficiency and numerical straightforwardness. Anisotropic geoelectrical structures are, however, often connected with spatially localized geological phenomena, such as shear zones, fractured complexes, etc., which do not conform the idea of a gradual, diffuse transition to the neighbouring units. Moreover, experiments with 2-D anisotropic models have shown that interactions of largely anisotropic domains with different anisotropy directions can result in strong distortions of the MT impedances (Pek and Verner, 1997). It is, therefore, desirable to limit both the spatial extent of the anisotropic domains and the anisotropy variations within those domains by employing modified regularization strategies which give preference to the ‘nonsmooth’ structural features. In this study, we have tested several regularization approaches, both quadratic and non-quadratic, which have recently been used in other geophysical interpretations (Portniaguine and Zhdanov, 1999). A summary of the regularization functionals used in our experiments is given in Table 1. We refer the reader to Portniaguine and Zhdanov (1999) for the mathematical background and the proof that the considered non-quadratic functionals are stabilizers in sense of the Tikhonov regularization theory. One of the serious questions of the regularized inversion (3.1) is the particular choice of the regularization parameter λ. One of the most complete analyses of this problem in a context of non-linear ill-posed geophysical inverse problems has been given by Haber and Oldenburg (2000), who proposed the generalized crossvalidation procedure as an efficient tool for selecting the global regularization parameter. Their inverse algorithm proceeds on an iteration by iteration basis, and the optimum regularization weight is modified in each linearized iteration step. We have not used this approach only because, at a later stage of the anisotropic inversion, we need to deal with more than only one penalty function, specifically two, one for the structure and another

144

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

Table 1 Summary of the regularization functionals used in our inversion tests for anisotropic layered conductivities Penalty

Theoretical formula



Numerical representation

∞



2

[m(z) − mprior (z)] dz

Least-squares norma

0 ∞ 

0

i

|∇m(z)| dz

(mi − mi−1 )2 + β2 , β → 0

0



∞

Minimum supportd

i

2

 Minimum gradient supporte

∞

[∇m(z)]2 2

0



2

[m(z) − mprior (z)]

[m(z) − mprior (z)] + β2

0

(mi − mi−1 )2



∞

Total variationc

)

i

[∇m(z)]2 dz

Roughnessb

prior 2

(mi − mi

[∇m(z)]

+ β2

dz, β → 0

i



dz, β → 0

i

prior 2

(mi − mi

prior 2 (mi − mi )

)

+ β2

(mi − mi−1 )2 (mi − mi−1 )2 + β2

,β→0

,β→0

a

Tikhonov and Arsenin (1977); minimizes the LS norm of the deviation of the model from a prior distribution. Constable et al. (1987); minimizes the roughness of the parameters in the LS sense. Prefers slow variations, tends to oversmooth the solution and suppress any sharp changes in the solution. c Rudin et al. (1992); minimizes the total variation of the parameter. It is indifferent to the size of the parameter jumps, and can preserve sharp edges in the solution. It is non-differentiable at points with zero gradient, hence introduction of a small stabilizing β is necessary. d Portniaguine and Zhdanov (1999); for a sufficiently small stabilizing β, it penalizes any deviation from mprior (z) by the same amount, independently of the size of the deviation. e Portniaguine and Zhdanov (1999); for a sufficiently small stabilizing β, it penalizes any sharp change in the solution, independently of the size of the jump. Tends to produce models with a minimum number of layers. b

one for the anisotropy, and the cross-validation becomes very computer intensive in this case. Another popular approach to choosing the regularization parameter is the L-curve criterion (Hansen, 1992). The L-curve is a plot of the norm of the regularized solutions Φstruct (m, mprior ) versus the norm of the residuals Φdat (m). This dependence has often an L-shaped form which reflects the heuristics that for large λ’s the residual increases without reducing the structural norm of the solution much, while for small λ’s the norm of the solutions increases rapidly without much decrease in the data residual. Thus, the best regularization parameter should lie on the corner of the L-curve. Though the L-curve approach is not flawless (see e.g., Engel and Grever, 1994), we have used it in our experiments to estimate the quasi-optimal regularization parameter λ. As in all non-linear inverse problems the stability of the inverse solution with respect to the regularization parameter is not guaranteed, the L-curve may become seriously malformed due to visiting different local minima of the target function when minimizing it with different λ’s. We further partly follow the approach by Kathirgamanathan et al. (2003) who suggest to start constructing the non-linear L-curve from a very small regularization parameter and increase it slowly, with always the preceding inverse solution taken as a starting guess for the next inversion run with the subsequent λ, so that

jumps between different modi of the target function (3.1) can be avoided as much as possible. 3.2. Local minimization Commonly, Newton-type methods are employed for minimizing the target function (3.1), particularly if quadratic regularization is implemented (e.g., Constable et al., 1987). The use of non-quadratic regularization functionals can, however, result in serious non-linearities of the target function and its derivatives, and may degrade, or even completely corrupt, the convergence of the Newton-type minimization methods. Portniaguine and Zhdanov (1999) propose to minimize (3.1) by using the conjugate gradient reweighted optimization with line search. They minimize the target function of a particular form

2



Φ(m, λ) = WM (m)(m − mprior )



2



+ λ−1 WD [dmod (m) − dobs ] → min, (3.2) where the specific structural norm given by the first term on the r.h.s. of (3.2) is shown to be compatible with all the regularization functions in Table 1. The minimization process in the iteration step (k + 1) is based on a

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

line search for the minimum of Φ[m(k) − τg(m(k) )] with respect to τ, where g(m(k) ) is the conjugate gradient direction, g(m(k) ) = g0 (m(k) ) + γk g(m(k−1) ), with γk



g0 (m(k) ) 2 =

,

g0 (m(k−1) ) 2 and g0 (m(k) ) is the gradient of (3.2) in the kth step, with the generally non-linear matrix WM (m) fixed with the arguments of the current iteration step, i.e., m = m(k) , g0 (m(k) ) = 2WTM (m(k) )WM (m(k) )(m(k) − mprior ) T

+ 2λ−1 [Smod (m(k) )] WTD WD [dmod (m(k) ) − dobs ], where Smod (m) =
M [dmod (m)]T is a matrix of the Frechet (parametric) derivatives of the forward solution. The structural weighting matrix WM (m) is updated in every iteration, which explains the reweighting character of the algorithm. The convergence of the conjugate gradient optimization has been found reliable in practical tests, with typically several tens of iteration steps required for the algorithm to converge to the solution for reasonable models. Another minimization procedure we have tested here is an analogy to the lagged diffusivity method that was proposed earlier for problems in image processing with the total variation regularization (Vogel and Oman, 1998). Let us express the parametric gradient of the structure norm in the form
M Φstruct (m) = 2WMD (m)(m − mprior ), which can be easily verified for all structure norms in Table 1. Then, the necessary condition for the target function (3.1) to reach a minimum is T

∇M Φ(m, λ) = 2λ−1 [Smod (m)] WTD WD [dmod (m)−dobs ] + 2WMD (m) (m − mprior ) = 0,

(3.3)

145

The convergence of the process (3.4) is generally not guaranteed for non-linear problems, and we have encountered situations in numerical tests when the parameter correction given by (3.4) resulted in increase of the target function between successive iterations. In such a case, the correction step was reduced by a predefined damping factor, repeatedly if necessary, until a decrease of the target function was achieved. With this modification, the algorithm converged reliably, with typically less than fifty iteration steps required to reach the solution. In Figs. 1 and 2, we demonstrate the performance of the inverse solution for a synthetic data set generated by a model with parameters given in Table 2. The model consists of five layers, two of which are anisotropic with different anisotropy strikes. The conductance of the deeper layer is chosen in such a way that the layer can be sensed in the surface data in the direction of its minimum resistivity. The model has been already analyzed earlier, specifically in the context of testing effects of a depth variable structural strike on MT data and results of MT decomposition procedures in Kov´acˇ ikov´a and Pek (2002). For the inversion tests, the impedance values were contamined with Gaussian noise with the standard deviation corresponding to 2% of the maximum impedance element for each period. The inverse solution was regularized by employing the standard roughness penalty applied to all three types of the model parameters involved, i.e., to the minimum and maximum resistivities and the anisotropy azimuth, with the same regularization weights used. This approach is clearly oversimplified as different penalization of the individual parameter types would be more appropriate, especially if the dynamic ranges of the anisotropic conductivity parameters differ significantly over the depths considered. Balancing several regularization weights would, however, increase the amount of computations significantly, so we have omitted the multiple weight regularization in our experiments to date. Fig. 1a shows the L-curve for various regularization parameters λ from 0.01 up to 1000. The values of λ in the knee region of the L-curve are between about 3 and 10. The plots in Fig. 1b show the resulting inverse models for λ = 0.316 (overfitting, model M1a), λ = 10 (close to the corner of the L-curve, model M1), and λ = 316.2 (oversmoothing, model M1b), respectively. For reference, the true model is also shown.

which is similar to g0 (m) above, except that now the matrix WMD (m) corresponds to the exact gradient of the respective structure norm. The lagged diffusivity method consists in applying the standard Gauss–Newton optimization procedure to (3.3) with the weighting matrix WMD kept fixed with the parameters of the current iteration step. Thus, we obtain a modified Gauss–Newton iteration process  −1 T m(k+1) = m(k) − [Smod (m(k) )] WTD WD Smod (m(k) ) + λWMD (m(k) )   T × [Smod (m(k) )] WTD WD [dmod (m(k) ) − dobs ] + λWMD (m(k) )(m(k) − mprior ) .

(3.4)

146

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

Fig. 1. (a) L-curve for the inversion of the synthetic data generated by the model given in Table 2, with the roughness penalty employed. The models corresponding to the points M1a, M1, and M1b are shown in the respective model panels in subfigure (b) along with the true model. (b) True model and anisotropic 1-D models from the L-curve in subfigure (a). In each pair of plots, the left plot shows the range between the minimum and maximum resistivities as a function of depth, the right plot shows the anisotropy azimuth.

Considering the model M1 in Fig. 1b as a solution to the regularized inverse problem, we can see that the two anisotropic conductors of the original synthetic model are correctly recovered as to their depth ranges (3–10 and 70–200 km), and quite satisfactorily also as regards the resistivity values (3 and 300 m for the shallow layer, and 30 and 300 m for the deep conductor). The shallow strike (−50◦ ) is recovered correctly as well, while the interpreted deep strike (about 30◦ ) is slightly greater than the true value (20◦ ). The most serious discrepancy between the original model and the inverse solution is the excessive anisotropy, of about one order of magnitude in terms of the max/min resistivity ratio, in the depth range between the two anisotropic conductors, which results from a poor resolution of the MT data with respect to the resistive layer (1000 m) within the range of 10–70 km. In subsequent simulations, we have tried to suppress the redundant anisotropy in the structure by imposing an additional constraint on the resistivity anisotropy within the model. The particular penalty function has

∞ been chosen Φanis = 0 log2 [ max (z)/ min (z)]dz (l2norm of the log-anisotropy) with a weighting factor wa . To find the model parameters m, we now minimized an extended target function Φ(m, λ, wa ) = Φstruct (m, mprior ) + wa Φanis (m) + λ−1 Φdat (m) → min .

(3.5)

In the first step, we have fixed the weight of the roughness penalty at λ = 10 and run the inversions with different values of the anisotropy penalty wa so that the product λwa ∈ (0.01, 1000). In this way, we could analyze how much of the anisotropy can be removed from the model M1 without affecting the data fit too much. The corresponding L-curve is displayed in Fig. 2a. The curve shows that for λwa greater than about 10, the residual norm increases rapidly, whereas for those values below 5 the fit is practically the same as for the original model obtained without any anisotropy penalty. The inverse solution for λwa = 5.62 is shown in Fig. 2b as a model M2.

Table 2 Parameters of a synthetic anisotropic five-layer model with two anisotropic layers with discordant anisotropy strikes Layer

h (km)

1 2 3 4 5

3 7 60 130 ∞

β S (◦ )

min (m)

−50

3

max (m)

Smin (S)

300

23.3

1000

3

1000 20

30

2333.3 60

300 200

Smax (S)

433.3

4333.3

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

147

Fig. 2. Top panels: (a) L-curve for a fixed roughness and a variable anisotropy penalty weight. The model corresponding to the point M2 is shown in subfigure (b). (c) Roughness vs. anisotropy dependence for models with fixed roughness penalty λ and variable anisotropy penalty wa . Bottom panels: (d) L-surface for variable roughness and anisotropy penalty weights. All (λ, λwa ) points used for constructing the contours are posted on the plot (empty circles). Three models from the knee region of the L-surface are shown as models M2a, M2, and M2b in Fig. 1b. The dashed contour marks the level of RMS = 1. The path corresponding to the L-curve from Fig. 1a (no anisotropy penalty, variable λ) is marked by gray large circles. The path corresponding to the L-curve from panel a (fixed λ = 10, variable wa ) is shown by large black square symbols. (e) Anisotropic 1-D models for various combinations of the roughness and anisotropy penalty weights. Positions of the models on the L-surface are indicated on the contour plot in panel (d).

A more objective approach is to employ both the structure and anisotropy penalty simultaneously via a multi-objective minimization procedure applied to the function (3.5). For this purpose, we plotted a 2-D constraint trade-off surface in Fig. 2d which generalizes the single-constraint L-curve to multiple constraints, and shows the data residuals as a function of both the roughness and anisotropy penalty weights. Though not L-shaped any more, the term L-surface, or hypersurface for more dimensions, has been in use for this constraint trade-off surface since Belge et al. (2002). The knee region of the L-surface is clearly delineated by a rapid change in the contours density and is close to the

isoline of Φdat (m) = 208 (shown by a bold dashed line in the plot), which corresponds to RMS = 1.0. This coincidence is not casual, and reflects the synthetic data generation process. Two models corresponding to the weights from the knee region of the L-surface are presented in Fig. 2e as models M2a and M2b, and they both show a relatively successful recovery of the anisotropic domains. As the periods used cover the section within the depth range of approximately 1.5–200 km, the anisotropies observed in the very shallow and very deep parts of the models have no significance. From comparing the above models for various regions of the knee zone of the L-surface, we may conclude that

148

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

the anisotropy penalty is constructive up to a certain limit only. If this limit is exceeded, the anisotropy decreases further, but at the cost of increasing the complexity of the structure. Fig. 2c demonstrates this dependency by showing the anisotropy versus roughness plot for our synthetic data set for a particular structural penalty weight λ = 10 and variable wa . In practice, models with an overestimated anisotropy penalty tend to produce suites of largely anisotropic individual peaks on an isotropic background (see model M2b in Fig. 2e), rather than giving clearly developed, compact anisotropic domains. In practice, models with more anisotropy, but with a simpler, smooth structure, will be likely better acceptable than models consisting of a series of thin anisotropic layers distributed throughout the model, though both may produce almost the same misfit. Then, our earlier approach to choosing the penalty weights individually one by one, i.e., first the roughness weight alone and only then the anisotropy penalty for the target function with the structure penalty factor fixed, may be applicable and computationally undoubtedly more economical. This approximate approach can be also qualitatively justified by observing the L-paths that correspond to the above two steps of the selection of λ and wa . In Fig. 2d, we additionally plotted the L-path for wa = 0 and a variable λ. From the corner of that path (model M1) we started another L-path, now with a fixed value of λ = 10 and a variable wa . The latter path closely follows the knee zone of the full L-surface, and only deviates significantly from the maximum curvature zone for very large anisotropy penalty weights, typically λwa > 30. The finally accepted model M2 is well within the knee zone of the L-surface.

Fig. 3 shows results of analogous experiments with the total variation penalty. The stabilizing parameter β was chosen 0.1 in this case. As an anisotropy penalty, the l1-norm of the log-anisotropy, Φanis = ∞ |log[ max (z)/ min (z)]| dz, was employed in these 0 inversion runs. The convergence properties of this approach are similar to those observed earlier with the roughness penalty used. However, the total variation minimization shows a better ‘de-smoothing’ performance as compared with the results in Fig. 2 and produces more focused image of the anisotropic layers. A more strongly non-linear structure penalty, particularly the minimum gradient support, has been also tested. Due to a highly irregular distribution and clustering of points that form the support for the corresponding Lsurface, optimal regularization weights are sometimes difficult to determine in this case. The minimum gradient support regularization is thus a tool more suitable for focusing and sharpening models that have been preinverted with some more conservative structural penalization, such as roughness, curvature or total variation, rather than for running the complete inversion from a fairly arbitrary initial model guess. 3.3. (Mis)Interpreting 2-D profile data by fitting 1-D anisotropic models Though exact indicators of the electrical anisotropy in the Earth’s subsurface are difficult to define, a quasiuniform split of perpendicular MT curves, particularly that of the MT phases, accompanied by negligible induction arrows across an area are often interpreted as due to the presence of deep anisotropic layers in

Fig. 3. L-surface and a representative inverse model for the synthetic test data from the model in Table 2. The dashed contour marks the level of RMS = 1. Large gray circles show the L-path for a variable λ and no anisotropy penalty applied (wa = 0). Large black squares mark the L-path for a fixed structure weight of λ = 5.62 and a variable anisotropy weight wa . A multi-objective regularization with a total variation penalty and the minimum l1-norm of the log-anisotropy is applied.

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

the Earth. Sometimes similar areal patterns of the MT functions may be, however, observed above isotropic but laterally inhomogeneous structures, especially if deep and large-scale heterogeneities are involved that do not produce significant induction arrows in a certain period range. In some cases, the vertical magnetic field may simply not be available. In this context, a natural question arises as to a possible misinterpretation of the electrical section if the split of the MT curves is falsely attributed to a deep electrically anisotropic structure.

149

To illustrate the character of possible errors that may result from formally interpreting a laterally non-uniform structure by a series of stitched 1-D anisotropic sections, we generated two synthetic data sets from two simple 2-D block models. The first model is an exact replica of a single-conductive-box-in-a-resistive-host-medium structure used in McNeice and Jones (2001) for testing their generalized MT decomposition scheme (see Fig. 4a for the model parameters). The other model is a semi-negative of the previous structure, and consists of a conductive layer interrupted by a resistive gap of a

Fig. 4. (a) Model of a conductive block (size 25 km × 4 km, depth 5 km, resistivity 50 m) in a resistive host medium, resistivity 1000 m, with marked sites at which the 1-D inversion was carried out. (b) Simulated experimental (diamonds) apparent resistivities and MT phases compared with curves obtained from the 1-D anisotropic inverse model (lines) at site G above the centre of the anomaly. Left panels: xy (full diamonds/full line) and xx (empty diamonds/dashed line) curves, right panels: yx (full diamonds/full line) and yy (empty diamonds/dashed line) curves. (c) Formal 1-D inverse anisotropic models at sites A through G. The top panels show ranges between the minimum and maximum resistivites, the bottom panels are for the anisotropy azimuths. The total variation regularization and the l2-norm log-anisotropy penalty were used, with weights typically λ ∈ (10.0, 17.8) and wa ∈ (1.0, 3.2).

150

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

Fig. 5. (a) Model of a conductive layer (thickness 4 km, depth 5 km, resistivity 50 m) interrupted by a resistive gap (width 25 km, resistivity 1000 m), with marked sites at which the 1-D inversion was carried out. Resistivity of the uniform host medium is 1000 m. (b) Simulated experimental (diamonds) apparent resistivities and MT phases compared with curves obtained from the 1-D anisotropic inverse model (lines) at site G above the centre of the anomaly. Left panels: xy (full diamonds/full line) and xx (empty diamonds/dashed line) curves, right panels: yx (full diamonds/full line) and yy (empty diamonds/dashed line) curves. (c) Formal 1-D inverse anisotropic models at sites A through G. The top panels show ranges between the minimum and maximum resistivites, the bottom panels are for the anisotropy azimuths. The total variation regularization and the l2-norm log-anisotropy penalty were used, with weights typically λ ∈ (10.0, 17.8) and wa ∈ (1.0, 3.2).

size identical to that of the conductive anomaly in the first model (Fig. 5a). The examples are rather artificial, as it would not be difficult to identify both the models as 2-D structures according to the spatial variations of their impedances and induction arrows, but we prefer to use them in this exercise in order to show the contrast between the true structures and their formal interpretations by 1-D anisotropic models more clearly. Before the interpretation was started, the synthetic impedances were rotated through −30◦ to produce nonzero secondary impedance elements, and then con-

tamined with Gaussian noise with the standard deviation corresponding to 2% of the maximum impedance element for each period. Diamond symbols in Figs. 4b and 5b show the simulated experimental apparent resistivities and phases for the site G situated in both cases on the surface above the centers of the anomalous domains. The 1-D inversions for anisotropic conductivities were carried out for a series of seven sites, A through G in Figs. 4a and 5a, covering a section of 32 km from the centers of the anomalous domains towards the left hand margin of the models. The period range of the simulated data

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

was 0.01–1000 s. We used the total variation regularization, as it seems to separate layers with largely different electrical parameters more sharply than the quadratic smoothing, and the l2-norm log-anisotropy penalty in the inverse runs. For the model with the conductive anomalous block, the interpreted 1-D model sections are displayed in Fig. 4c. For the sites A through C, they first show a gradually emerging, slightly anisotropic layer with the strike close to 30◦ . Starting from the site D, the first one situated above the anomaly, the real conductive anomaly is recovered relatively well, both as to its position and resistivity. It shows only casual anisotropy that does not seem to be significant. Above the anomaly, however, a false, largely anisotropic block is interpreted, with the strike of −60◦ , reaching down to the bottom of the model. This anisotropic structure accounts for the split of the MT curves and cannot be reduced nor eliminated from the 1-D models. The fit of the 1-D model data to the experimental impedances deteriorates towards the center of the anomaly, from 1.11 at site A to 3.00 at G, in terms of the RMS. Nonetheless, visually, the coincidence between the experimental and model data seems to be quite satisfactory even for the worst case of site G (lines in Fig. 4b), and the most serious discrepancies are found for secondary impedances in a range of periods where their experimental modules are very small. For the model of the resistive gap within a conductive layer (Fig. 5), we observe similar spurious anisotropy effects in formal 1-D interpretations. Above the normal part of the section, between sites A and C, the conductive layer is recovered relatively well, with a slightly anisotropic shadow zone beneath the true position of the conductor. Above the resistive gap, the MT curves split can be only explained by a largely anisotropic two-level structure, with roughly the ratio of 300 m/1000 m and azimuth of 30◦ within the depth range of 8–30 km, and 1000 m/10,000 m, azimuth 30◦ , at greater depths. Again, this large anisotropy cannot be eliminated from the models. Similarly as above, the fit between the experimental and 1-D model data deteriorates when the center of the anomaly is approached, starting from RMS = 1.59 for site A and ending with RMS = 2.34 at G. To give a chance to inspect the misfit visually, the 1-D model curves for the worst case of site G are again plotted in Fig. 5b along with the simulated experimental data. Both the above examples demonstrate that modelling of 2-D MT curves by stitched 1-D anisotropic models may result in serious misinterpretation of the electrical structure of the underground. Though the detailed features of the artifact structures may be different depending on the particular 2-D data generating model, attempts

151

to account for the split of the MT curves due to 2-D non-uniformity by using 1-D models always result in relatively large spurious anisotropies, especially in the deep part of the section. These effects may be obviously potentially dangerous when interpreting data over large-scale heterogenous structures that are not sufficiently covered by measurements and do not provide a clear pattern of the induction arrows. Though an elemetary exercise, the above analysis may become relevant if tests for true anisotropic structures in laterally inhomogeneous settings are performed, as, e.g., with the practical MT data from the southern Portugal in Section 3.4 below. As interactions of lateral conductivity inhomogeneities with discordant anisotropic structures may result in serious distortions of MT fields (e.g., Pek and Verner, 1997), assessing information obtained from stitched 1-D anisotropic images of such complex models may become important in a situation when a true inverse procedure for laterally inhomogeneous conductivity distributions with arbitrary anisotropy is still not available. In the ensuing experimental analysis, we dealt with specific questions of (i) how much the apparent anisotropy due to lateral conductivity variations changes if the deep anomaly is overlain by a shallower anisotropic crustal screen, as well as (ii) what artifacts can be expected to occur when the properties of the anisotropic screen are estimated from local 1-D models. The synthetic tests we performed were based on the two 2-D models, a single anomalous conductive box and a resistive window, presented above. In both cases, the lateral anomaly was covered by an anisotropic layer, 3 km thick, with the resistivities 10 m/100 m and with the anisotropy strike of 30◦ with respect to the axis of homogeneity of the 2-D structure. A third model was also considered, identical to the conductive block anomaly but with the resistivity of the block decreased by a factor of 10, i.e., to 5 m. The generated synthetic MT curves were again rotated by −30◦ and contamined with 2% of Gaussian noise to preserve conditions comparable to those of the previous experiments. In Fig. 6, we present 1-D inverse anisotropic models for the above 2-D structures at the central point above the anomalous domain (corresponding to G in Figs. 4 and 5). For all three cases, the presence of the anisotropic layer results in reduced apparent anisotropy due to the lateral inhomogeneities. It is clearly an effect of the increased conductivity of the shallow screening layer and does not originate from its anisotropy. As to its position, depth, max/min resistivities and strike, the anisotropic layer is recovered correctly in all three cases. A conductive anomaly beneath the anisotropic layer is resolved satisfactorily as to its minimum resistivity and its position,

152

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

Fig. 6. Formal 1-D anisotropic inversions of synthetic MT curves from the site G above the center of the 2-D models from Figs. 4 and 5 with the inhomogeneous layer covered by an anisotropic layer (anis) with the thickness of 3 km, min/max resistivities of 10 m/100 m and anisotropy strike of 30◦ with respect to the direction of homogeneity. The total variation regularization and the l2-norm log-anisotropy penalty were used in the inversion. (a) Results for the model with a conductive block with the resistivity of 50 m (see Fig. 4) covered by the anisotropic layer. Parameters of the inversion: λ = 56.2, wa = 0.56, RMS = 1.76. (b) As model (a), only the resistivity of the anomalous block is 5 m. Parameters of the inversion: λ = 100.0, wa = 1.0, RMS = 4.30. (c) Results for the model with a resistive gap in a conductive layer of 50 m (see Fig. 5) covered by the anisotropic layer. Parameters of the inversion: λ = 17.8, wa = 1.76, RMS = 1.44.

with 1–2 km error in its depth estimate at the most. Even in the case of the 50 m block (Fig. 6a), with a total conductance of 80 S, the resistivity estimate between 25 and 35 m is quite good, considering that the conductance of the anisotropic cover is 300 S/30 S in the directions of the max/min conductivity, respectively. The maximum resistivity within the conductive anomalous block is less well constrained by the inversion, and the 1-D resistivity image may suggest that only one continuous anisotropic layer with depth-variable anisotropy occupies the whole depth range between 2 and 10 km. The underlying resistive substratum is recovered relatively sharply and shows only moderate apparent anisotropy aligned with the structural strike. From the above synthetic experiments we can conclude that 1-D images of shallow anisotropic layers are not severely affected by lateral conductivity anomalies situated at greater depths and can be recovered well by a 1-D inversion. Due to the effect of the anisotropic layer, however, the resolution with respect to the deeper con-

ductivity anomalies is reduced in a direction-dependent way (compare with the effect of the polarized spectacles by Jones, 2006), which can further result in apparently merging both the anisotropic layer and the underlying anomaly into one single anisotropic structure, with not a clearly marked interface between the vertical substructures. Of course, this relatively simple picture may become obscured if near-surface distortions, both static and inductive, affect the MT data to a considerable degree. 3.4. Inversion of experimental MT curves The developed inverse algorithm has been applied to field data from MT experiments in southern Portugal, in the transition of two southernmost prominent zones of the Iberian Variscides, the Southern Portuguese Zone (SPZ) and the Ossa Morena Zone (OMZ). Recent studies have shown that a long and narrow belt of mafic and ultramafic rocks lines up at the border between these

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

zones. Those rocks are currently interpreted as an ophiolite complex (Badajoz Acebuches Ophiolite Complex, BAOC) placed along a suture formed during accretion of the exotic terranes of the SPZ onto the autochthonous or parautochthonous terranes (Ribeiro et al., 1990; Quesada et al., 1994; Fonseca and Ribeiro, 1993). In the Portuguese mainland the BAOC represents the southern border of the Beja Igneous Complex (BIC) deformed in the vicinity of the Ferreira-Ficalho thrust (Figueiras et al., 2002). The Ferreira-Ficalho thrust (FF) is a left-handed, south-thrusting WNW-ESE tectonic structure that runs continuously, except when interrupted by the Ficalho fault, between the Tagus and Guadalquivir Cainozoic basins, separating the SPZ and OMZ (Figueiras et al., 2002). It is known that a local increase of the shear stress and/or a notable change in fluid chemistry, by means of mixing processes, when associated with considerable volume of aqueous-carbonic fluids can produce high graphite concentrations (e.g., Large et al., 1994). Micro and macro-channels developed in shear zones are characterised by preferential directions and, consequently, the electrical conductivity distribution in these zones may exhibit anisotropy. The geological context of the study area is then greatly favourable to the presence of anisotropic structures at upper to middle crustal depths. References to the presence of anisotropic black schists in south Iberia, as well as to the anisotropic behaviour of MT data can be found in Pous et al. (2004). Independent data from seismic refraction profiles carried out in the SPZ show evidence of seismic velocity anisotropy (Matias, 1997). Matias’ (1997) interpretation of the seismic data in terms of an anisotropic variation gives a fast direction of N125◦ E, or N55◦ W, with the maximum velocity value of 6.4 km/s, probably due to foliation in schistose formations. An extensive broad-band geoelectrical induction experiment in the south of Portugal was carried out along a more than 200 km long profile crossing the main Variscan zones in the region. A detailed report on the results can be found in Almeida et al. (2001). Here we only concentrate on a short section of the profile, with five MT sites over a distance of about 25 km, situated above the old suture zone between the SPZ and OMZ. The original 2-D REBOCC inversion (Siripunvaraporn and Egbert, 2000) of the MT curves from the whole profile produced a model with a very good fit to the experimental data (Almeida et al., 2001). In the BAOC/BIC region above the suture zone, the conductivity model is marked by three prominent electrically anomalous structures. Specifically, they are (i) a shallow conductor in the south (labeled C2 in Almeida et

153

al., 2001) that correlates with zones generated by the Variscan collision, BAOC and IOMZOS (Internal Ossa Morena Zone Ophiolitic Sequences), (ii) an extensive conductive body in the middle to lower crust (C3), correlated with the old suture between the SPZ and OMZ, and (iii) an upper crustal resistor (R1), correlated with the BIC. We show the approximate position of those structures by dashed lines and corresponding labels later in Fig. 8. Though an excellent fit was achieved by 2-D inverting the principal MT curves, some specific features in the data suggest that non-negligible non-2-D effects influence the MT field. Specifically, variations of single-site, single-frequency regional strike estimates (Fig. 7) for the sites considered show rather a systematic frequency dependence. While they are poorly defined and practically random for short periods up to about 0.1 s, they all show a similar systematic trend for longer periods, traversing a broad range of directions, exceeding even 60◦ at some sites (09, 01, 03). At the two northernmost sites, 09 and 01, penetration of a discordant structure may be indicated by strike excursions at periods between 1 and 100 s. Though a reasonable regional strike estimate of about N100◦ E, or N80◦ W, has been suggested as a result of a multi-site, multi-frequency decomposition for the whole profile in (Almeida et al., 2001), the systematic variations of the strike at our selected sites indicate a locally more complex structure with possibly anisotropic or discordant domains. Fig. 7 shows 1-D anisotropic inverse models for the five selected sites. To constrain the structure and anisotropy, we used the roughness and the l2-norm of anisotropy penalties in the target function (3.5). The particular parameters of the inversion runs are given in the caption to Fig. 7 along with the RMS of the final models. The RMS misfits are in general rather large in these inversions, which is mainly due to a poorer fit to the experimental secondary impedances. As experimental data rarely meet exactly the basic 1-D anisotropy condition, Zxx = −Zyy (see, e.g., Abramovici, 1974; Kov´acˇ ikov´a and Pek, 2002), the inversion necessarily produces models with secondary impedances fitting some average of the real secondary curves. Extreme cases are highly distorted curves for which one of the principal impedance phases leaves its standard quadrant in some range of directions, as is the case with the curve 03 in our data set. The 1-D models show a rather consistent anisotropic domain starting at shallow depths of 1–2 km at the northern sites and less than 1 km further to the south. The dominant anisotropy is observed at upper crustal depths, between about 3 km down to about 10–20 km. Considering the results of the synthetic modelling in the previous

154

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

Fig. 7. Top panel: single site, single period Monte Carlo estimates of the regional strike at five sites along a S–N profile in the transition zone between the SPZ and OMZ in southern Portugal. The size of the bubbles is inversely proportional to the RMS for the resulting composite models. Bottom panel: formal 1-D anisotropic inversions of the MT data from the same five sites. The roughness penalty and l2-norm of log-anisotropy penalty were used in the inversion. Parameters of the inversion were for all cases λ ∈ (56.2,100.0) and wa ∈ (0.56, 1.0). RMS misfits for the inverse models are 7.49, 17.02, 3.01, 4.62, 2.19 for the sites 04 through 09, respectively.

section, this anisotropic domain may well accomodate the large conductor C3 from the original 2-D REBOCC model without requiring its anisotropy (see Fig. 8). Considering only the dominant features of the inverse 1-D sections, we furher constructed a simplified 2-D model for the selected MT sites by stitching the 1-D models together. The aim of the 2-D anisotropic modelling is to explain the systematic non-2-D features in the MT data in the area of study. Since this aim is still

far from being completed, we only present in Fig. 8 one tentative model from a suite of experiments based on a trial-and-error modification of the starting stitched structure. In the middle panel of the figure, we show the fit of the principal resistivity and phase curves to the experimental data. The fit is still rather qualitative for the present model, and especially discrepancies in the secondary impedances (not shown) are still unacceptably large. Nonetheless, there are some positive aspects even

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

155

Fig. 8. Top panel: example of a tentative 2-D model across the old suture zone between the SPZ and OMZ in southern Portugal, obtained by trial-anderror adjusting the stitched 1-D anisotropic sections from Fig. 7. The hatched zones indicate anisotropic domains with the anisotropy strikes defined with respect to the structural strike of the model, which was chosen N70◦ W in this particular region. The bold dashed lines approximately delineate anomalous resistivity zones C2, C3, and R1 from the REBOCC inversion by Almeida et al., 2001. Middle panel: fit between the principal apparent resistivity and phase curves for the above model and experiment for the sites 04 through 09. The curves are shown in the NS–EW coordinates. Gray circles and full lines: NE experimental and model curves, respectively; empty circles and dashed lines: EN experimental and model curves, respectively. Bottom panel: comparison of single site, single period Monte Carlo regional strike estimates from the experimental impedances (empty circles) and from the data produced by the model in the top panel (gray circles) for the sites 04 through 09. The model data were contamined with 2% of Gaussian noise for the stochastic sampling. The size of the bubbles is inversely proportional to the RMS for the resulting composite models, the reference bubble is shown beneath the plots. The strikes are related to the geographic north.

156

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

in the present version of the model that can be summarized as follows: 1. As compared with the 1-D models, the anisotropic domain could be reduced to only the upper crustal part of the 2-D model, with an isotropic conductor at greater crustal depths, which can be correlated with the conductor C3 of Almeida et al. (2001). The position of the interface between the two domains is not resolved unambiguously. 2. The anisotropy strike of 15◦ (with respect to the structural strike of the model, i.e., about N55◦ W in geographical coordinates) is in good agreement with the crustal seismic anisotropy reported by Matias (1997). 3. The model can account for the excessive split between the NS and EW curves at sites 03 and 04 in the south of the mini-profile. As a part of this split results from the energy transfer from the B-mode (electric currents perpendicular to the structural strike) into the E-mode (currents parallel to the structural strike) via intermode coupling due to the anisotropy, it can be controlled by relatively small variations of the anisotropy parameters of the crustal layer. 4. Even for the realistic anisotropy strike of only 15◦ , the model can produce relatively large and frequency dependent variations in the apparent regional strike estimates on the surface, showing strike excursions as large as 60◦ , though not yet identical with those observed for the experimental impedances. To illustrate the frequency dynamics of the regional strike estimates for the present model structure, we have plotted the experimental strikes jointly with the corresponding estimates for the model structure in the bottom panel of Fig. 8. The model strikes have been estimated from impedances contamined artificially with 2% Gaussian noise, which roughly corresponds to the experimental data variances at periods greater than 1 s. Since we still have not attempted to model the shallowest part of the structure in detail, the short period strike estimates indicate a well-developed 2-D structure that nearly copies the strike of the crustal anisotropic layer. For periods greater than about 1 s, the composite model does not seem to be a satisfactory approximation to the model impedances. The formal strike estimates show large directional variations for these periods, similar to those observed in the experimental data, especially at sites 02 and 03. For periods longer than about 100 s, the experimental and model strike estimates diverge, and further modelling is still necessary to explain the turn of the experimental strikes in this period range.

4. Conclusion MT inversion for anisotropic conductivities bears some specific features as compared to isotropic models. In this paper, we have presented a complete algorithm for a 1-D MT inversion for anisotropic layered conductors, which represents a straightforward extension of the standard Occam idea (Constable et al., 1987) to anisotropic media. Though more advanced 2-D, and even 3-D inverse codes have been newly under study in geoelectrical applications, e.g. in Pain et al. (2003), many specific features of the inversion for anisotropic conductivities can be still better understood and analyzed in a simple 1-D context. As MT data alone cannot principally resolve the complete conductivity tensor in laterally uniform structures, the inverse procedure has been developed for a simplified 1-D anisotropic model with only the azimuthal anisotropy considered within the layers. Additional information has to be drawn from non-MT data in order to assess further details of the anisotropic conductivity structure, as, e.g., the resolution of the anisotropy dip (e.g., Yin, 2003). The solution to the direct problem and the calculation of parametric sensitivities is based on a stable; impedance propagation method for 1-D anisotropic conductors adopted from Pek and Santos (2002). For the inversion, a straightforward extension of the Occam algorithm by Constable et al. (1987) to the anisotropic case is used, with the parameter set modified to involve the three independent anisotropy conductivity parameters in each layer, specifically the minimum and maximum resistivities and the anisotropy azimuth. As the anisotropy in the Earth is often connected with spatially sharply demarcated geological structures, we tested various regularization approaches within the inverse algorithm, with both quadratic and non-smooth structure complexity penalties employed. Introducing an additional penalty on the overall resistivity anisotropy shows that a large mass of the anisotropy predicted by unconstrained inversion can be removed, or largely reduced, without affecting the model-data fit significantly. Thus, efficient focusing on anisotropic domains can be achieved by the bi-criterial penalization. The minimization of the non-linear target function was based on either the conjugate gradient reweighted optimization (Portniaguine and Zhdanov, 1999), or on the analogy to the lagged diffusivity iterative solution (Rudin et al., 1992). Both methods perform well with nonsmooth regularizators, the former being more reliable as to the convergence, but slower than the lagged diffusivity method.

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

We used a rather qualitative version of the L-curve criterion (Hansen, 1992) for estimating quasi-optimal weights for the regularization penalties in our 1-D anisotropic inverse approach. A proper selection of the regularization parameters can be partly objectivized by using that tool, though the choice of the regularization weights still depends on subjective views, especially in the version with the multi-criterion regularization for both the structure complexity and anisotropy. Experience with using the developed inverse algorithm showed its stable and reliable performance for synthetic models with multiple discordant anisotropic domains situated at different depths. Problems with properly selecting the regularization weights can sometimes arise if non-quadratic regularization penalties in the target functional (3.5) are used. In those cases, the points that form the support of the L-surface may show a non-regular distribution and/or clustering, which makes it difficult to reconstruct the topology of the L-surface. Then, a two-step approach to the inversion is recommended, with the quadratic regularization applied in the initial stage of the minimization, and with subsequent parameter focusing through the use of a non-quadratic complexity penalization. Using 1-D anisotropic inversion in the modelling of field MT data is only exceptionally justifiable by a truly layered anisotropic setting. A truly 1-D anisotropic Earth requires that the induction parameters show a sufficiently uniform distribution across an area and that no significant vertical magnetic fields are induced. A good example of a successful interpretation by a 1-D anisotropic model has been recently presented by Linde and Pedersen (2004) for a radio-MT survey over limestone formations. In the majority of practical situations, however, the apparent anisotropy of the MT curves represents either a manifestation of lateral inhomogeneities within the Earth, or of their combination with anisotropic conductors, or an effect of static distortions. Then, the 1-D anisotropic inversion can be only justified as a simple imaging tool, which gives the possibility to process simultaneously all the impedance information provided by the complete impedance tensor. The risk is that largely anisotropic spurious structures may be indicated in the 1-D models, as demonstrated by our synthetic experiments in Section 3.3. We have observed that stitched 1-D anisotropic images of a 2-D conductive or resistive box can recover the true structure satisfactorily down to the depth of the true anomalous domain. Beneath the anomalies, however, a substratum with excessive spurious anisotropy is interpreted in both the conductive and resistive case. As a true 2-D MT inversion for arbitrarily anisotropic conductivities is still not available, employing the 1-D

157

anisotropic imaging is useful in designing 2-D, or even 3-D, test models for anisotropy simulation studies. As a practical example, the crust close to an old suture zone between the South Portuguese Zone and Ossa Morena Zone in southern Portugal bears both geological and geophysical indications on possible anisotropy due to shearing and graphitization (Santos et al., 1999; Almeida et al., 2001). For a small subset of MT data from this region, we suggest an extended interpretation consisting of a combination of a strike analysis, 1-D anisotropic inversions and a trial-and-error 2-D anisotropic modelling. Though still at a very preliminary and qualitative stage, the results obtained promise an improvement to explaining some systematic, truly non-2-D features in the experimental data, like, e.g., systematic trends in the frequency dependence of regional strike estimates, or phase excursions out of their natural quadrants. Acknowledgements The authors wish to thank Randall Mackie and two anonymous reviewers for their helpful comments and criticism on the earlier versions of the manuscript, as well as to the Guest Editor Alan Jones for his help. The financial support provided by the Grant Agency of the Academy of Sciences of the Czech Republic, through the grant No. A3012401, and by the Czech Sci. Found., through the grant No. 205/04/0746, is gratefully acknowledged. The authors would like to thank the Academy of Sciences of the Czech Republic and the GRICES Portugal for promoting the Czech-Portuguese cooperation within the project ’Inversion of electromagnetic depth sounding data for heterogeneous and anisotropic conductivities in the Earth’. References Abramovici, F., 1974. The forward magnetotelluric problem for an inhomogeneous and anisotropic structure. Geophysics 39, 56–68. Abramovici, F., Shoham, Y., 1977. Inversion of anisotropic magnetotelluric data. Geophys. J. R. Astr. Soc. 50, 55–74. Almeida, E., Pous, J., Monteiro Santos, F., Fonseca, P., Marcuello, A., Queralt, P., Nolasco, R., Mendes-Victor, L., 2001. Electromagnetic imaging of a transpressional tectonics in SW Iberia. Geophys. Res. Lett. 28, 439–442. Belge, M., Kilmer, M.E., Miller, E.L., 2002. Efficient determination of multiple regularization parameters in a generalized L-curve framework. Inverse Problems 18, 1161–1183. Constable, S.C., Parker, R.L., Constable, C.G., 1987. Occam’s inversion: a practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics 52, 289–300. Dekker, D.L., Hastie, L.M., 1980. Magnetotelluric impedances of an anisotropic layered Earth model. Geophys. J. R. Astr. Soc. 61, 11–20.

158

J. Pek, F.A.M. Santos / Physics of the Earth and Planetary Interiors 158 (2006) 139–158

Eisel, M., Haak, V., 1999. Macro-anisotropy of the electrical conductivity of the crust: a magnetotelluric study from the German Continental Deep Drilling site (KTB). Geophys. J. Int. 136, 109–122. Engel, H.W., Grever, W., 1994. Using the L-curve for determining optimal regularization parameters. Numer. Math. 69, 25–31. Figueiras, J., Mateus, A., Gonc¸alves, M.A., Waerenbourgh, J., Fonseca, P., 2002. Geodynamic evolution of the South Variscan Iberian Suture as recorded by mineral transformations. Geodinamica Acta 15, 45–61. Fonseca, P., Ribeiro, A., 1993. Tectonics of the Beja-Acebuches ophiolite: a major suture in the Iberia Variscan foldbelt. Geol. Rundsch. 82, 440–447. Haber, E., Oldenburg, D., 2000. A GCV based method for nonlinear ill-posed problems. Comput. Geosci. 4, 41–63. Hansen, P.C., 1992. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev. 34, 561–580. Jones, A.G., 2006. Electromagnetic interrogation of the anisotropic Earth: Looking into the Earth with polarized spectacles. Phys. Earth Planet. Inter. 158, 281–291. Kathirgamanathan, P., McKibbin, R., McLachlan, R.I., 2003. Source release rate estimation of atmospheric pollution from a non-steady point source, Part 2: Source at an unknown location. Res. Lett. Inf. Math. Sci. 5, 85–118. Kov´acˇ ikov´a, S., Pek, J., 2002. Generalized Riccati equations for 1-D magnetotelluric impedances over anisotropic conductors, Part I: Plane wave field model. Earth Planets Space 54, 473–482. Large, D.J., Christy, A.G., Fallick, A.E., 1994. Poorly crystalline carbonaceous matter in high grade metasediments: implications for graphitisation and metamorphic fluid compositions. Contrib. Mineral. Petrol. 116, 108–116. Li, X.B., Oskooi, B., Pedersen, L.B., 2000. Inversion of controlledsource tensor magnetotelluric data for a layered Earth with azimuthal anisotropy. Geophysics 65, 452–464. Linde, N., Pedersen, L.B., 2004. Evidence of electrical anisotropy in limestone formations using the RMT technique. Geophysics 69, 909–916. Loewenthal, D., Landisman, M., 1973. Theory for magnetotelluric observations on the surface of a layered anisotropic halfspace. Geophys. J. R. Astr. Soc. 35, 195–214. Mareschal, M., Kellett, R.L., Kurtz, R.D., Ludden, J.N., Bailey, R.C., 1995. Archean cratonic roots, mantle shear zones, and deep electrical anisotropy. Nature 373, 134–137. Matias, L., 1997. The Contribution of Experimental Seismology to the Knowledge of the Crustal Structure in Portugal Mainland. Ph.D. thesis, University of Lisbon (in Portuguese). McNeice, G.W., Jones, A.G., 2001. Multisite, multifrequency tensor decomposition of magnetotelluric data. Geophysics 66, 158–173. O’Brien, D.P., Morrison, H.F., 1967. Electromagnetic fields in an Nlayered anisotropic half-space. Geophysics 32, 668–677. Pain, C., Herwanger, J.V., Saunders, J.H., Worthington, M.H., Oliveora, C.R.E., 2003. Anisotropic resistivity inversion. Inverse Problems 19, 1081–1111.

Pek, J., Verner, T., 1997. Finite-difference modelling of magnetotelluric fields in two-dimensional anisotropic media. Geophys. J. Int. 128, 505–521. Pek, J., Santos, F.A.M., 2002. Magnetotelluric impedances and parametric sensitivites for 1-D generally anisotropic layered media. Comput. Geosci. 28, 939–950. Pous, J., Munoz, G., Heise, W., Melgarejo, J.C., Quesada, C., 2004. Electromagnetic imaging of Variscan crustal structures in SW Iberia: the role of interconneted graphite. Earth Planet. Sci. Lett. 217, 435–450. Portniaguine, O., Zhdanov, M.S., 1999. Focusing geophysical inversion images. Geophysics 64, 874–887. Quesada, C., Fonseca, P., Munh´a, J., Oliveira, J., Ribeiro, A., 1994. The Beja-Acebuches Ophio-lite (southern Iberia Variscan Foldbelt): geological characterization and geodynamic significance. Boletin Geol´ogico y Minero de Espana 105 (1), 3–49. Reddy, I.K., Rankin, D., 1971. Magnetotelluric effect of dipping anisotropies. Geophys. Prosp. 19, 84–97. Regis, C.R.T., Rijo, L., 1997. 1-D inversion of anisotropic magnetotelluric data. Extended Abstracts Book, 50th IC Soc. Bras. Geof. Ed., Brazil, vol. II, pp. 673–674. Regis, C.R.T., Rijo, L., 2000. Approximate equality constraints in the inversion of anisotropic MT data. Abstracts Book, 15th Workshop on Electromagnetic Induction in the Earth, Cabo Frio, Brazil, p. 47. Ribeiro, A., Quesada, C., Dallmeyer, R.D., 1990. Geodynamic evolution of the Iberian Massif. In: Dallmeyer, R.D., Martinez Garcia, E. (Eds.), Pre-Mesozoic Geology of Iberia. Springer-Verlag, pp. 398–409. Rudin, L.L., Osher, S., Fatemi, E., 1992. Nonlinear total variation based noise removal algorithms. Physica D 60, 259–268. Santos, F.A.M., Mendes-Victor, L.A., 1997. 1-D interpretation of anisotropic MT data. Extended Abstracts Book, 50th IC Soc. Bras. Geof. Ed., Brazil, vol. II, pp. 772–775. Santos, F.A.M., Pous, J., Almeida, E.P., Queralt, P., Marcuello, A., Matias, H., Mendes-Victor, L.A., 1999. Magnetotelluric survey of the electrical conductivity of the crust across the Ossa Morena Zone and South Portuguese Zone suture. Tectonophysics 313, 449–462. Siripunvaraporn, W., Egbert, G., 2000. An efficient data-subspace inversion method for 2D magnetotelluric data. Geophysics 65, 791–803. Tikhonov, A.N., Arsenin, V.Y., 1977. Solution of ill-posed problems. W.H. Winston and Sons, Washington D.C, 272 pp. Vogel, C.R., Oman, M.E., 1998. Fast, robust total variation-based reconstruction of noisy, blurred images. IEEE Trans. Image Proc. 7, 813–824. Wannamaker, P.E., 2005. Anisotropy versus heterogeneity in continental solid Earth electromagnetic studies: fundamental response characteristics and implications for physicochemical state. Surv. Geophys. 26, 733–765. Yin, C.C., 2003. Inherent nonuniqueness in magnetotelluric inversion for 1D anisotropic models. Geophysics 68, 138–146.