**Multi-site, multi-frequency tensor decomposition of magnetotelluric data**

*Gary W. McNeice*, Phoenix Geophysics Ltd.
Alan G. Jones, Geological Survey of Canada*

**Summary**

Accurate interpretation of magnetotelluric data relies upon an
understanding of the directionality and dimensionality inherent
in the data, and correct implementation of a method for
removing the effects of small-scale galvanic scatterers. The
galvanic distortion analysis advocated by Groom and Bailey
has become the *de facto* standard. and rightly so given that the
approach both decomposes the magnetotelluric impedance
tensor into determinable and indeterminable parts, and tests
statistically the validity of the approximation. In its original
formulation, one fits the model on a frequency-by-frequency
and site-by-site basis independently, and determines
qualitatively the most consistent regional strike angle. Whilst
this approach has the attraction that one gains a more intimate
understanding of the dataset, it is time-consuming and requires
repetitive, tedious application. We propose an extension to
Groom-Bailey decomposition in which a global minimum is
sought for the most appropriate strike direction and distortion
parameters for a range of frequencies and a set of sites. We
illustrate application of the analysis to synthetic datasets.

**Introduction**

In recent years, impedance tensor decomposition analysis has
become an integral part of modern interpretation of
magnetotelluric (MT) data. The objective of this analysis is to
remove the first-order effects of galvanic distortion produced
by near surface zones of anomalous conductivity, and to
recover the impedances and strike of the regional geoelectric
structure(s). Of the parameters derived during decomposition
analysis, the regional strike is often the poorest resolved
parameter (Groom *et al.* 1993; Jones & Groom 1993), which
presents somewhat of a problem, since the accuracy of the
interpretation is dependent on the correct adoption of a regional
strike consistent with the data. In this GB analysis, galvanic
distortion is characterized by two decomposition parameters
*shear* and *twist*, and induction by the *regional strike*. These
three are iteratively constrained to find frequency-independent
estimates required by the decomposition model (see Groom *et
al.* 1993, for a tutorial). Regional strikes determined at each
frequency of a set of sites are then compared to obtain an
estimate of regional strike for the entire data set (assuming the
data are two-dimensional and sensing the same regional 2D
structure). For a more representative estimate of regional
strike, one really needs to perform a weighted average of the
determined strike estimates, since the resolution of regional
strike is a function of frequency, site location, and data
precision.

In this paper, an extension to Groom-Bailey decomposition is proposed in which the physical constraints implicit to the decomposition model are imposed simultaneously on the data set and a global minimum sought. This procedure enables the determination of the estimate of strike most consistent with the data set in a statistical sense. The variation of strike resolution as a function of frequency and position are used in the determination of the strike estimate. The approach additionally provides the interpreter with information on the degree to which the data fits the model of regional two dimensionality with galvanic distortions.

Synthetic and real data are used to demonstrate the robustness
of the proposed decomposition analysis.

**Groom-Bailey Decomposition**

The basic model is of non-inductive, galvanic distortion by small-scale three-dimensional (3D) bodies in thin surface sheet over a two-dimensional (2D) regional Earth (3D/2D)

(Bahr 1984), where **R** is a rotation tensor. The elements of the
2x2 tensor **C** are real and frequency-independent. There have
been many approaches proposed to derive the parameters
descriptive of the regional conductivity structure; which are the
geoelectric strike, , of **R**, and the two complex impedances of
**Z _{2D}.** The most appealing to date is the approach of Groom &
Bailey (1989), who describe a tensor decomposition based on
factoring the telluric distortion tensor

where **g** is a scalar termed "site gain" and **T**, **S**, and **A** are
tensor factors termed *twist, shear *and *anisotropy *respectively.

The twist, **T,** and shear, **S,** tensors are the determinable parts of
the distortion matrix, and the site gain **g** and anisotropy tensor
**A** form the indeterminable parts.

The seven parameters of the GB factorization are the two
descriptors of telluric distortion twist **t** and shear **e**, the regional
strike , and the four parameters contained in the two complex
regional impedances, Z_{xy} and Z_{yx}. The parameters of the
factorization are found by statistically fitting the seven
parameters of the decomposition model to the eight (8) data of
the measured impedance tensor.

In formal GB analysis the decomposition model is fit to the
measured data frequency by frequency and site by site.
However, for the 3D/2D distortion model to be valid frequency
independent estimates of the twist, shear and regional strike
must be found over a band of frequencies and a group of
neighbouring sites. In practice this is done by iteratively
constraining the estimates of twist, shear and strike found at
each frequency in an attempt to find a frequency band for
which the constrained values result in an acceptable misfit error
(Groom *et al.* 1993). If the regional conductivity structure is
truly 1D or 2D, the model will hold in general at low
frequencies reducing the task to finding an upper frequency
limit below which the distortion becomes solely galvanic
(distorter becomes inductively small) and the distortion
parameters become frequency independent.

**Regional strike estimation**

Of the parameters solved for in decomposition analysis of the impedance tensor, regional strike is by far the most important parameter. Accurate recovery of the regional impedances, and therefore accuracy of the interpretation, is dependent on recovery of the correct regional strike. Unfortunately, regional strike is often one of the poorest resolved parameters in decomposition analysis (Jones & Groom 1993). The resolution of regional strike is a function of both frequency and position and its recovery is impeded by galvanic distortion and experimental error.

The importance of accurate recovery of the regional strike can be seen by examining properties of the measured impedance tensor under Groom-Bailey factorization. In an arbitrary coordinate system, the measured impedance tensor is given by (1). The elements of the impedance tensor are distortion and rotation dependent mixtures of the two regional impedances. In the regional coordinate system however, (1) reduces to

and the two columns of the impedance tensor contain scaled
estimates of the regional impedances (denoted by the prime).
The impedance elements contained in the two columns are
linearly dependent, having the same impedance phase, differing
solely in magnitude. This phase relationship was first
recognized by Bahr (1984). The magnitude of the four
elements is described by the twist (**t**) and shear (**e**) factors of
the GB factorization.

Equation (3) shows that theoretically scaled estimates of the regional impedances can be found in the presence of galvanic distortion by simply rotating the impedance tensor to the regional coordinate system. An error in the regional strike estimate will results in the recovery of impedances which are mixtures of the regional impedances, while errors in twist and shear result solely in errors in the magnitudes of the recovered regional impedances. Recovery of the correct regional strike is therefore essential for accurate recovery of the regional impedances.

In the decomposition procedures proposed by Zhang *et al.*
(1987), Bahr (1988), Chakridi *et al.* (1992), and Smith (1995)
the linear dependence of the columns of the impedance tensor
is used to find the regional strike. Their procedures solve for
a strike angle for which the elements of the two columns have
equal impedance phases. This method however is unstable in
the presence of noise since depending on distortion (twist and
shear) some of the elements of the tensor can be small in
magnitude and dominated by error, producing significant errors
in the recovered regional strike (Jones & Groom 1993).

Groom-Bailey decomposition performs better in the presence
of noise since it simultaneously solves for all of the parameters
of the decomposition model. Twist and shear are often better
resolved than the regional strike as they are dependent on the
relative magnitude of the impedance elements as opposed to
phase properties of the elements. The inclusion of twist and
shear into a model fitting all eight elements of the impedance
tensor stabilizes the estimation of regional strike.

**Extension of Groom-Bailey Decomposition for Multiple
Sites and Multiple Frequencies**

To improve the estimation of regional strike, one can extend the GB analysis to fit statistically an entire dataset simultaneously with the assumed 3D/2D distortion model. The regional strike is usually only well defined over a small subset of the data., and thus this extension yields superior strike estimation since it accounts for the position and frequency dependence of strike resolution. The dependence of model misfit on strike defines a regional strike statistically consistent with the entire data set.

In the extended GB analysis, **S** magnetotelluric sites having **N**
frequencies (**S*N*8** data, for data at every frequency from
every site) can be fit with **S*(N*4 + 2) + 1** unknowns. The
unknowns are the **S*N** regional complex impedances (Z_{xy} and
Z_{yx}), the **S*2** descriptors of telluric distortion (twist **t** and shear
**e**) and the regional strike angle (). In the extended model the
**S** sites are required to have the same regional strike; the telluric
distortion parameters shear and twist remain site dependent
although independent of frequency. Assuming that all the data
are independent, the model has **S*(N*4 - 2) - 1** degrees of
freedom. In the case of a single site, the extended model is
equivalent to that of the GB method with the strike, twist and
shear constrained to be independent of frequency, although the
solutions are arrived at differently.

The model parameters are found in the extended analysis by
minimizing a misfit function based on the observed and model
summary decomposition coefficients. The elements of each of
the modelled impedance tensors are described by four complex
nonlinear equations, and the minimization procedure seeks a
simultaneous solution to the **S*N*8** nonlinear equations (4 real
and 4 imaginary equations for each of the **S*N** impedance
tensors) describing the measured data set. An objective
functional normalized by the data variances was chosen in
order to reduce any bias produced by frequency and position
dependence of experimental error and the frequency
dependence of impedance magnitudes. Additionally, the
objective function allows one to derive significance levels for
the parameter estimates obtained by the minimization
procedure.

Errors on the parameters are derived using a bootstrap method
similar to that of Groom & Bailey (1991).

**Failure of decomposition analysis**

There are two situations that can occur in decomposition
analysis of the impedance tensor, where the method fails to
recover the regional impedances accurately. Decomposition
fails when the shear approaches 45^{0} (**e** of 1) or when the sum
of distortion and regional anisotropy approaches an **s** of 1. In
both situations, decomposition analysis fails because the
measured impedance tensor approaches singularity and the
decomposition model becomes underdetermined.

In the case of high shear the two rows of the impedance tensor become linearly dependent in any coordinate system. The strike therefore is unresolvable, as the impedance tensor assumes the form of a distorted 2D tensor in any coordinate system. The shear angle will be well resolved by the linear dependence of the two rows, although only the sum of the regional strike and twist angles will be resolved. In order to recover the regional impedances the strike must be found from adjacent sites or other independent information.

In the high anisotropy case, the distortion parameters cannot be
resolved and while a strike angle is resolved, it cannot be
determined if it is that of the regional structure or a local
distorter. If the strike determined by decomposition is in
agreement with the regional strike found at other nearby sites
(in an inductive sense) it may be reasonable to assume that the
recovered impedances are those of the regional structure.

**Examples of Application**

**Synthetic Example I**

To illustrate the performance of the proposed decomposition
procedure in the presence of experimental noise, thirty-one
realizations of the theoretical impedance tensor discussed by
Jones & Groom (1993) were generated. The distorted tensor
has a site gain (g) of 1.06, an anisotropy (s) of 0.172, a twist
angle of -2.1^{} (t = -0.037), and a shear angle of 24.95^{} (e =
0.47). The inductive 2D regional strike is 0^{}.

Gaussian noise, having a standard deviation of 4.5% of the magnitude of the largest impedance element, was added to the tensor (4) to produce thirty-one noise-contaminated realizations. These thirty-one realizations form an ideal data set to test the performance of the decomposition procedure in the presence of noise, as all the realizations are equally sensitive to the inductive strike and obey the telluric distortion model differing solely by experimental noise.

Figure 1 shows histograms of the strike, twist and shear found from an unconstrained GB analysis of the thirty-one realizations. The distributions of the GB parameters clear shows that shear and twist are more stable under decomposition than strike angle determination. The strike directions determined appear to be approximately normally distributed, although they ill define the regional strike.

The results of the extended decomposition analysis yields a
twist angle of -2.2^{}, a shear angle of 25.6^{} and strike direction
of 0.1^{}, very close to the correct values of -2.1^{}, 24.95^{} and 0^{}
respectively. The results appear to be accurate and free of
noise bias, in contrast to the outcome of iteratively constraining
the frequency-dependent decomposition parameters as
advocated by Groom et al. (1993). Such a procedure resulted
in an estimate of strike of -7^{} 4^{} (Jones & Groom 1993).

**Synthetic Example II**

Now we evaluate decomposition analysis of a more realistic theoretical response, where sensitivity to the 2D inductive strike is a function of frequency. Groom & Bailey (1991) produced an accurate 3D/2D dataset by superposing 2D numerical and 3D analytic responses. The 3D analytical response is that of a small conducting hemisphere, and the 2D host is a fault-like structure. For the frequencies considered the hemisphere has a negligible inductive response, but it has a significant effect on both the electric and magnetic fields of the regional 2D structure. For frequencies <10 Hz (periods >0.1 s), the distortion effects produced by the anomalous magnetic fields of the hemisphere become negligible and the hemisphere acts as a galvanic scatterer.

The strike of the regional 2D structure is 30^{}. Groom *et al.*
(1993) have shown that the strike is poorly defined at both
short and long periods. For short periods the data are
insensitive to the 2D boundaries of the regional structure and
are essentially 1D, whereas at long periods the 2D structure
becomes inductively thin, having only a galvanic response
which combines with the galvanic response of the
hemispherical distorter. The phases of the two regional
impedances become virtually identical, and the regional
response becomes an anisotropic 1D response (2D galvanic
distortion of the response of the deep regional 1D structure).

To simulate experimental noise, 2% Gaussian noise was added
to the theoretical responses. Decomposition analysis, prior to
the addition of noise, yields a twist angle of -12^{} and shear
angle of 30^{} (Groom *et al.* 1993). The 3D scattering effects of
the hemisphere in this data set are first-order effects while the
inductive response of the 2D structure has a much smaller
effect, making this a good dataset for testing the extended
decomposition routine.

The results of the extended decomposition analysis are shown
superimposed on the results of the unconstrained GB analysis
in Fig 2. The extended analysis yields a strike direction of
30.7^{}, twist angle of -13.0^{} and shear of 30.6^{}. The extended
model has 69 degrees of freedom, with a 95% confidence level
having a ** ^{2}** value of 90.0. The chi-squared misfit of the model
is 49.5, well below the 95% confidence level.

The extended analysis recovers each of the regional responses
to within a multiplicative static shift factor. The magnitude and
phase of the Z_{yx} (TM mode) impedance are contaminated by
noise at long periods where the Z_{yx} impedance is small in
magnitude. This noise contamination is enhanced by the
distortion anisotropy, which reduces the magnitude of the Z_{yx}
impedance. However, the long period Z_{yx} impedances would
have the highest scatter even in the absence of distortion due to
their small magnitude.

**Conclusions**

We have extended the Groom-Bailey decomposition method to
find the most consistent 2D parameters from a set of sites over
a given frequency band. This method allows one to determine
rapidly a regional strike statistically consistent with frequency,
site location and data precision dependencies. The recovered
impedances, where the 3D/2D decomposition model holds, are
free of first-order effects of galvanic distortion, and the
procedure provides a validation for further 2D interpretation.
However, we advocate careful use of this global minimization,
as the parameters may vary with frequency as, for example, the
strike may be depth-dependent (see, e.g., Marquis *et al.*, 1995).

**References**

Bahr, K. 1984. Elimination of local 3D distortion of the
magnetotelluric tensor impedance allowing for two different
phases. Contributed paper at *Seventh Workshop on
Electromagnetic Induction in the Earth and Moon*, held in
Ile-Ife, Nigeria, August 15-22.

Bahr, K. 1988. Interpretation of the magnetotelluric impedance
tensor, regional induction and local telluric distortion. *J.
Geophys.*, **62**, 119-127.

Bahr, K. 1991. Geologic noise in magnetotelluric data: a
classification of distortion type. *Phys. Earth. Planet. Inter.*,
66, 24-38.

Chakridi, R., Chouteau, M. & Mareschal, M. 1992. A simple
technique for analysing and partly removing galvanic distortion from the magnetotelluric impedance tensor: application to Abitibi and Kapuskasing data (Canada).
*Geophys. J. Inter.*, **108**, 917-929.

Groom, R.W. & Bailey, R.C. 1989. Decomposition of
magnettotelluric impedance tensors in the presence of local
three-dimensional galvanic distortion. *J. Geophys. Res.*, **94**,
1913-1925.

Groom, R.W. & Bailey, R.C. 1991. Analytic investigations of
the effects of near-surface three-dimensional galvanic
scatters on MT tensor decompositions. *Geophysics*, **56**, 496-518.

Groom, R.W., Kurtz, R.D., Jones, A.G. & Boerner, D.E., 1993.
A quantitative methodology for determining the
dimensionality of conductivity structure and the extraction
of regional impedance responses from magnetotelluric data.
*Geophys. J. Inter,* **115**, 1095-1118.

Jones, A.G. & Groom, R.G. 1993. Strike angle determination
from the magnetotelluric impedance tensor in the presence
of noise and local distortion: rotate at your peril! *Geophys.
J. Inter.*, **113**, 524-534.

Marquis, G., Jones, A.G. & Hyndman, R.D., 1995. Coincident
conductive and reflective lower crust across a thermal
boundary in southern British Columbia, Canada. *Geophys.
J. Int.*, **20**, 111-131.

Smith, J.T., 1995, Understanding telluric distortion matrices.
*Geophys. J. Int.*, **122**, 219-226.

Zhang, P., Roberts, R.G. & Pedersen, L.B. 1987.
Magnetotelluric strike rules. *Geophysics*, **52**, 267-278.