Appraisal of the Magnetotelluric and Magnetovariational Transfer Functions’ Selection in a 3-D Inversion

Abstract

Magnetotelluric (MT) and magnetovariational (MV) sounding are two principal geophysical methods used to determine the electrical structure of the earth using natural electromagnetic signals. The complex relationship between the alternating electromagnetic fields can be defined by transfer functions, and their proper selection is crucial in a 3-D inversion. A synthetic case was studied to assess the capacity of these transfer functions to recover the electrical resistivity distribution of the subsurface and to evaluate the advantages and disadvantages of using the tipper vector W to complement the impedance tensor Z and the phase tensor Φ. The analysis started with two sensitivity tests to appraise the sensitivity of each type of transfer function, which is calculated for an oblique conductor model, showing that the resistivity perturbation of the same model will produce distinct perturbations to different transfer functions; the transfer function sensitivity is significantly different. A 3-D inversion utilizing the quasi-Newton method based on the L-BFGS formula was performed to invert different transfer functions and their combinations, along with quantifying their accuracy. The synthetic case study illustrates that a 3-D inversion of either the Z or Φ responses presents a superior ability to recover the subsurface electrical resistivity; joint inversions of the Z or Φ responses with the W responses possess superior imaging of the horizontal continuity of the conductive block. The appraisal of the 3-D inversion results of different transfer functions can facilitate assessing the advantages of different transfer functions and acquiring a more reasonable interpretation.

1. Introduction

Two principal forms of geophysical methods used to determine the electrical structure of the earth through the natural existing electromagnetic fields have been in wide use since the last century, which consist of the magnetotelluric (MT) sounding method and the magnetovariational (MV) sounding method [1]. During an MT and/or MV investigation, the presence of a common model with a plane electromagnetic wave vertically incident on a plane horizontally inhomogeneous earth [2] allows one to simultaneously measure time variations in the natural electromagnetic field E_h (E_x, E_y) and H (H_x, H_y, H_z) at the earth’s surface [3]. The complex relationship between the alternating electromagnetic fields can be defined by transfer functions based on the Maxwell equations [4,5]. In principle, the basic transfer functions in the MT consist of the impedance tensor Z and the apparent resistivities ρ_ij calculated from the components Z_ij of the impedance tensor Z, whereas the basic transfer function in the MV is the tipper vector W (the Wiese–Parkinson vector). For completeness, we also review the interstation impedance tensor Q, which is defined from the relationship between the survey electric field and the reference magnetic field [6,7], as well as the interstation horizontal magnetic tensor M and the interstation horizontal electric (telluric) tensor T [8,9]. To recover the regional phase information directly from the observed impedance tensor, researchers also focus on investigating the phase tensor Φ [10], which is defined from the real and imaginary parts of the MT impedance tensor Z. We offer a remark concerning notations: we employ for the tipper vector the notation W following the original terminology [3] and use for the horizontal electric (telluric) tensor the notation T, which is always reserved for the tipper vector [11,12]. Since the current field data acquisition is still dominated by conventional MT and/or MV surveys, only the traditional MT and MV transfer functions are considered in this manuscript, including the Z, W, and Φ responses.

In the traditional scheme of the MT and MV data interpretation, estimation of the transfer functions at multiple survey sites and frequencies, with subsequent inversion of the estimated transfer functions in terms of the subsurface conductivity, constitutes the essence of this scheme [13,14]. With advancing computational resources of high-performance clusters, three-dimensional (3-D) inversion techniques have gained a lot of attention in recent years and are now a more widely used tool for the scheme mentioned above [15,16,17,18,19]. Since the wavelet representation of a signal is inherently multiscale [20,21,22], the sparsity constraint in the wavelet domain has also been utilized in the electromagnetic inversion to recover the multiscale features of the subsurface [23,24]. Compared to 1-D or 2-D inversion, the numbers of elements and degrees of freedom in a 3-D inversion are much larger and must be constrained by data [25], so that the decision regarding the transfer functions’ data selection is one important preparation in 3-D inversion. Therefore, it is necessary to compare the characteristics and sensitivities of each data type for the reasonable selection of different data types to obtain more realistic 3-D inversion results.

In the current 3-D inversion of MT data, as the on-diagonal components of the MT impedance tensor Z are nonzero in the 3-D environments and are strongly affected by off-profile geoelectrical structures (e.g., data are acquired predominantly along the 2-D profiles) to at least some degree [26,27], the approach of inverting the full components of the MT impedance tensor Z are hence used instead of inverting only the off-diagonal components (Zo, and the subscript “o” denotes the off-diagonal) in the 3-D inversion [28,29]. However, the use of on-diagonal components is sometimes worthless because of the quality reduction in the on-diagonal components of Z caused by the presence of anthropogenic noise [30,31]. In addition, there has been a widespread recognition that MT impedance tensor data are often distorted by unresolvable small-scale structures [32], on which the charges on conductivity boundaries or gradients can cause galvanic distortion and then generate electric field effects. However, the MV transfer functions W and M are unaffected by this galvanic distortion [33]. Hence, the approach of inverting the tipper vector W or the interstation horizontal magnetic tensor M individually and including W and M simultaneously when performing a 3-D inversion of impedance tensor Z have long been recognized as the available schemes to improve the MT data interpretation and constrain the subsurface [7,28,33,34,35]. Analogous to the basic MV transfer function (W), the phase tensor Φ is also unaffected by the galvanic distortion, even though it is derived from the observed impedance tensor Z [10,36]; moreover, it is a statically undistorted 3-D response [32]. Thus, the approach of inverting the phase tensor Φ for galvanically distorted data sets using 3-D inversion has been implemented in several publications [37,38,39,40].

As mentioned above, although there are many comparative studies on the inversion of different MT and MV data, there are few systematic analyses and research studies in view of the sensitivities of different data types. Therefore, we attempt to use multiple approaches to demonstrate the sensitivities of different data types, including the identification of the locations and periods responsive to the geoelectrical anomaly, the quantification of the anomalous effects on different data responses, and the calculation of different sensitivity functions. Based on this, we examine the advantages and disadvantages of a 3-D inversion of the traditional MT and MV transfer functions. This paper aims to appraise the ability of these transfer functions to recover the electrical resistivity distribution of the subsurface, as well as evaluating the advantages and disadvantages of combining the Z, Zo, and Φ responses with the W responses when modeling the characterization of the subsurface in a 3-D environment.

2. Methodology

2.1. The 3-D Electromagnetic Inversion Method

Among the most widely used methods for solving the optimization problem (e.g., Occam’s inversion method, Siripunvaraporn et al. [27]; the nonlinear conjugate gradients method, Newman and Alumbaugh [41]; the Gauss–Newton method, Sasaki [42]; the limited-memory quasi-Newton method, Newman and Boggs [43]), a limited memory quasi-Newton method based on the Broyden–Fletcher–Goldfarb–Shanno (BFGS) formula is preferred in this study. This method constitutes an approximation of the inverse of the Hessian matrix by defining the curvature pair, which does not need additional storage requirements, and thus effectively accelerates the convergence of the iterative inversion process. In the L-BGFS inversion, the product of “${\hat{\mathrm{H}}}_k$” and ${\mathrm{G}}_k$ is used to solve for the change in the model at iteration step k + 1, where “${\hat{\mathrm{H}}}_k$” is the approximation of the inverse of the Hessian matrix, which is updated at every iteration using the limited-memory BFGS formula [44], and ${\mathrm{G}}_k$ is the gradient vector of the objective functional. The implementation of the L-BFGS inversion is completed by using the ModEM software package [15,16] and is discussed in detail in Yu et al. [45].

2.2. The Explicit Forms of the Matrix L in Different Cases of Transfer Functions

The matrix L was introduced by Egbert and Kelbert [15] and represents the linearized data functional. It depends on the details of the observation functionals (e.g., the MT and MV transfer functions), as it is applied to the discrete numerical forward solution. In this manuscript, we attempt to give the explicit forms of the L operators in different cases of transfer functions and follow the notation and approach of Egbert and Kelbert [15]; the expression of the matrix L is expressed as

$${\mathit{L}}_j={\left.\frac{\mathrm{\partial }{\psi }_j}{\mathrm{\partial }\mathbf{e}}\right|}_{{\mathbf{e}}_{0,}{\mathbf{m}}_0}=\sum _{j=1}^k\frac{\partial {\psi }_j}{\partial e_j}.$$

(1)

As indicated in Equation (1), the expression of the matrix L is not the same for different nonlinear data functionals ${\psi }_j$. Consequently, to obtain the discrete numerical forward solution to different data functionals requires implementation of the derivation of the L operator arising in different transfer functions.

Under the assumption that the plane electromagnetic wave source is assumed to be two independent field sources that correspond to spatially uniform sources of a fixed frequency polarized in the orthogonal directions (the x- and y-directions), the MT impedance tensor Z defined from relations between the horizontal components of the electric and magnetic fields of point observation can be expressed as

$$\left(\genfrac{}{}{0pt}{}{{\mathbf{E}}_{\mathit{h}}^{\mathrm{SX}}}{{\mathbf{E}}_{\mathit{h}}^{\mathrm{SY}}}\right)=\left[\begin{array}{cc}{\mathrm{Z}}_{xx}& {\mathrm{Z}}_{xy}\\ {\mathrm{Z}}_{yx}& {\mathrm{Z}}_{yy}\end{array}\right]\left(\genfrac{}{}{0pt}{}{{\mathbf{H}}_{\mathit{h}}^{\mathrm{SX}}}{{\mathbf{H}}_{\mathit{h}}^{\mathrm{SY}}}\right),$$

(2)

where the superscripts “SX” and “SY” denote the field sources polarized in the x- and y-directions, respectively. From Equation (2), the rows of the impedance tensor Z take the following form

$$\begin{gathered} {\mathrm{Z}}_{xx}=\frac{{\mathrm{E}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{E}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}, \\ {\mathrm{Z}}_{xy}=\frac{{\mathrm{E}}_x^{\mathrm{SY}}{\mathrm{H}}_x^{\mathrm{SX}}-{\mathrm{E}}_x^{\mathrm{SX}}{\mathrm{H}}_x^{\mathrm{SY}}}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}, \\ {\mathrm{Z}}_{yx}=\frac{{\mathrm{E}}_y^{\mathrm{S}\mathrm{X}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{E}}_y^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}, \\ {\mathrm{Z}}_{yy}=\frac{{\mathrm{E}}_y^{\mathrm{SY}}{\mathrm{H}}_x^{\mathrm{SX}}-{\mathrm{E}}_y^{\mathrm{SX}}{\mathrm{H}}_x^{\mathrm{SY}}}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}. \end{gathered}$$

(3)

Together with Equation (1), the operator L_Zij can be explicitly expressed as

$$\begin{gathered} L_{{\mathrm{Z}}_{\mathrm{xx}}}=\frac{L_{e,x}({\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_y^{\mathrm{SX}})-L_{b,x}\left({\mathrm{E}}_x^{\mathrm{SX}}{-\mathrm{E}}_x^{\mathrm{SY}}\right)}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}{-\mathrm{Z}}_{xx}\frac{{L_{b,x}(\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_y^{\mathrm{SX}})-L_{b,y}({\mathrm{H}}_x^{\mathrm{SY}}{-\mathrm{H}}_x^{\mathrm{SX}})}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}, \\ {L}_{{\mathrm{Z}}_{xy}}=\frac{L_{e,x}\left({\mathrm{H}}_x^{\mathrm{SX}}-{\mathrm{H}}_x^{\mathrm{SY}}\right)-L_{b,x}\left({\mathrm{E}}_x^{\mathrm{SX}}{-\mathrm{E}}_x^{\mathrm{SY}}\right)}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}-{\mathrm{Z}}_{xy}\frac{{L_{b,x}(\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_y^{\mathrm{SX}})-L_{b,y}({\mathrm{H}}_x^{\mathrm{SY}}{-\mathrm{H}}_x^{\mathrm{SX}})}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}, \\ {L}_{{\mathrm{Z}}_{yx}}=\frac{L_{e,y}({\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_y^{\mathrm{SX}})-L_{b,y}\left({\mathrm{E}}_x^{\mathrm{SY}}{-\mathrm{E}}_x^{\mathrm{SX}}\right)}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}-{\mathrm{Z}}_{yx}\frac{{L_{b,x}(\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_y^{\mathrm{SX}})-L_{b,y}({\mathrm{H}}_x^{\mathrm{SY}}{-\mathrm{H}}_x^{\mathrm{SX}})}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}, \\ {L}_{{\mathrm{Z}}_{yy}}=\frac{L_{e,y}({\mathrm{H}}_x^{\mathrm{SX}}-{\mathrm{H}}_x^{\mathrm{SY}})-L_{b,x}\left({\mathrm{E}}_y^{\mathrm{SY}}{-\mathrm{E}}_y^{\mathrm{SX}}\right)}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}-{\mathrm{Z}}_{yy}\frac{{L_{b,x}(\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_y^{\mathrm{SX}})-L_{b,y}({\mathrm{H}}_x^{\mathrm{SY}}{-\mathrm{H}}_x^{\mathrm{SX}})}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}, \end{gathered}$$

(4)

where L_e denotes the L operator of the electric field component, and L_b denotes the L operator of the magnetic field component. The explicit form of the operator L_Zij given here facilitates efficiently calculating the Jacobian matrix.

Other basic transfer functions in the MT are the apparent resistivity and impedance phases, which are calculated from the components of the impedance tensor Z. Separating the complex impedance tensor into their real (X) and imaginary (Y) parts, the general expressions for the apparent resistivity and impedance phases can be written as

$${\rho }_{ij}=\frac{{\left|{\mathrm{Z}}_{ij}\right|}^2}{\omega {\mu }_0}=\frac{{\left(X_{ij}\right)}^2+{\left(Y_{ij}\right)}^2}{\omega {\mu }_0},$$

(5)

$${\phi }_{ij}={\tan }^{-1}\left(Y_{ij}/X_{ij}\right).$$

(6)

Since the apparent resistivity and impedance phases are real, whereas the impedance tensor is complex, we take the derivative of ρ_ij (Equation (5)) concerning the solution vector u using the chain rule, and we take only the real part. Then, we have

$$\begin{gathered} \frac{\partial {\rho }_{ij}}{\partial \mathbf{u}}=\frac{\partial {\rho }_{ij}}{\partial X_{ij}}\frac{\partial X_{ij}}{\partial \mathbf{u}}+\frac{\partial {\rho }_{ij}}{\partial Y_{ij}}\frac{\partial Y_{ij}}{\partial \mathbf{u}}=\frac{2}{\omega {\mu }_0}(X_{ij}\frac{\partial X_{ij}}{\partial \mathbf{u}}+Y_{ij}\frac{\partial Y_{ij}}{\partial \mathbf{u}}) \\ =\frac{2}{\omega {\mu }_0}[X_{ij}\mathrm{Re}(\frac{\partial {\mathrm{Z}}_{ij}}{\partial \mathbf{u}})+Y_{ij}\mathrm{Im}(\frac{\partial {\mathrm{Z}}_{ij}}{\partial \mathbf{u}})]=\mathrm{Re}(\frac{2{\mathrm{Z}}_{ij}^{\ast }}{\omega {\mu }_0}\frac{\partial {\mathrm{Z}}_{ij}}{\partial \mathbf{u}}) \\ =\mathrm{Re}(\frac{2{\mathrm{Z}}_{ij}^{\ast }}{\omega {\mu }_0}{\mathbf{L}}_{{\mathrm{Z}}_{ij}}). \end{gathered}$$

(7)

Hence, the matrices L of the apparent resistivity and impedance phase can be expressed as

$${\mathbf{L}}_{{\rho }_{ij}}=\frac{2{\mathrm{Z}}_{ij}^{\ast }}{\omega {\mu }_0}{\mathbf{L}}_{{\mathrm{Z}}_{ij}}{, \mathbf{L}}_{{\phi }_{ij}}=\frac{i{\mathrm{Z}}_{ij}^{\ast }}{{\left|{\mathrm{Z}}_{ij}\right|}^2}{\mathbf{L}}_{{\mathrm{Z}}_{ij}},$$

(8)

where the superscript asterisk (*) denotes the complex conjugate.

The phase tensor Φ is defined from the real and imaginary parts of the MT impedance tensor Z, which can be written as the matrix (Booker 2014)

$$\mathbf{\Phi }={\mathbf{X}}^{-1}\mathbf{Y}=\left[\begin{array}{cc}{\mathrm{\Phi }}_{xx}& {\mathrm{\Phi }}_{xy}\\ {\mathrm{\Phi }}_{yx}& {\mathrm{\Phi }}_{yy}\end{array}\right]=\frac{1}{\det \left(\mathbf{X}\right)}\left(\begin{array}{cc}{X}_{yy}{Y}_{xx}-X_{xy}{Y}_{yx}& X_{yy}{Y}_{xy}-X_{xy}{Y}_{yy}\\ X_{xx}{Y}_{yx}-X_{yx}{Y}_{xx}& X_{xx}{Y}_{yy}-X_{yx}{Y}_{xy}\end{array}\right),$$

(9)

where X and Y are the real and imaginary parts of the MT impedance tensor Z, respectively, and det(X) = X_xxX_yy− X_yxX_xy is the determinant of X. Given that the phase tensor is real, we deal with the derivative of the phase tensor Φ (Equation (9)), concerning the solution vector u using the chain rule, and we take only the real part

$$\begin{gathered} \frac{\partial \mathbf{\Phi }}{\partial \mathbf{u}}=\frac{\partial {(\mathbf{X}}^{-1}\mathbf{Y})}{\partial \mathbf{u}}=-{\mathbf{X}}^{-1}\frac{\partial \mathbf{X}}{\partial \mathbf{u}}{\mathbf{X}}^{-1}\mathbf{Y}+{\mathbf{X}}^{-1}\frac{\partial \mathbf{Y}}{\partial \mathbf{u}} \\ =-{\mathbf{X}}^{-1}{\mathbf{X}}^{-1}\mathbf{Y}\mathrm{Re}(\frac{\partial \mathbf{Z}}{\partial \mathbf{u}})+{\mathbf{X}}^{-1}{\mathbf{X}}^{-1}\mathbf{X}\mathrm{Im}(\frac{\partial \mathbf{Z}}{\partial \mathbf{u}}) \\ ={\left({\mathbf{X}}^{-1}\right)}^2(-\mathbf{Y}\mathrm{Re}(\frac{\partial \mathbf{Z}}{\partial \mathbf{u}})+\mathbf{X}\mathrm{Im}(\frac{\partial \mathbf{Z}}{\partial \mathbf{u}}\left)\right) \\ =\mathrm{Re}(-i{\left(\frac{\mathrm{adj}(\mathbf{X})}{\det (\mathbf{X})}\right)}^2{\mathbf{Z}}^{\ast }\frac{\partial \mathbf{Z}}{\partial \mathbf{u}})=\mathrm{Re}(-i{\left(\frac{\mathrm{adj}(\mathbf{X})}{\det (\mathbf{X})}\right)}^2{\mathbf{Z}}^{\ast }{\mathbf{L}}_{\mathrm{Z}}), \end{gathered}$$

(10)

where adj(X) is the adjugate of X. The matrices L of the phase tensor can be expressed as

$${\mathbf{L}}_{{\mathrm{\Phi }}_{ij}}=-i{\left(\frac{\mathrm{adj}(\mathbf{X})}{\det (\mathbf{X})}\right)}^2{\mathbf{Z}}^{\ast }{\mathbf{L}}_{{\mathrm{Z}}_{ij}}.$$

(11)

The basic transfer function in the MV is the tipper vector W; similar to the impedance tensor, it is defined from the relations between the vertical and horizontal magnetic fields. Then, the tipper vector W can be expressed as

$$\left(\begin{array}{cc}{\mathbf{H}}_z^{\mathrm{SX}}& {\mathbf{H}}_z^{\mathrm{SY}}\end{array}\right)=\left[\begin{array}{cc}{\mathrm{W}}_{zx}& {\mathrm{W}}_{zy}\end{array}\right]\left(\begin{array}{cc}{\mathbf{H}}_{\mathit{h}}^{\mathrm{SX}}& {\mathbf{H}}_{\mathit{h}}^{\mathrm{SY}}\end{array}\right),$$

(12)

where the rows of the tipper vector W take the following form

$$\begin{gathered} {\mathrm{W}}_{zx}=\frac{{\mathrm{H}}_z^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_z^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}, \\ {\mathrm{W}}_{zy}=\frac{{\mathrm{H}}_z^{\mathrm{SY}}{\mathrm{H}}_x^{\mathrm{SX}}-{\mathrm{H}}_z^{\mathrm{SX}}{\mathrm{H}}_x^{\mathrm{SY}}}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}. \end{gathered}$$

(13)

Then, similar to the operator L_Zij, the operator L_Wij can be explicitly expressed as

$$\begin{gathered} L_{{\mathrm{W}}_{zx}}=\frac{L_{b,z}\left({\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_y^{\mathrm{SX}}\right)-{\mathrm{L}}_{b,y}\left({\mathrm{H}}_z^{\mathrm{SY}}-{\mathrm{H}}_z^{\mathrm{SX}}\right)}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}-{\mathrm{W}}_{zx}\frac{{L_{b,x}(\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_y^{\mathrm{SX}})-L_{b,y}({\mathrm{H}}_x^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SX}})}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}, \\ {L}_{{\mathrm{W}}_{zy}}=\frac{L_{b,z}({\mathrm{H}}_x^{\mathrm{SX}}-{\mathrm{H}}_x^{\mathrm{SY}})-L_{b,x}({\mathrm{H}}_z^{\mathrm{SX}}-{\mathrm{H}}_z^{\mathrm{SY}})}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}-{\mathrm{W}}_{zy}\frac{{L_{b,x}(\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_y^{\mathrm{SX}})-{\mathrm{L}}_{b,y}({\mathrm{H}}_x^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SX}})}{{\mathrm{H}}_x^{\mathrm{SX}}{\mathrm{H}}_y^{\mathrm{SY}}-{\mathrm{H}}_x^{\mathrm{SY}}{\mathrm{H}}_y^{\mathrm{SX}}}. \end{gathered}$$

(14)

3. Synthetic Model Study

To conduct the synthetic model study, we first calculated the forward responses from a synthetic model with an embedded oblique conductor and then generated a sensitivity test to evaluate the effects of the oblique conductor on the different types of transfer functions data by using the calculated responses, which consisted of Zo, ρ, Z, W, and Φ. The sensitivity function calculations were also performed, unraveling the sensitivities for different responses. Finally, we performed the 3-D inversion based on the L-BFGS algorithm of different transfer functions data individually; moreover, we implemented the 3-D joint inversion with different combinations of Zo, Z, or Φ with W, which aimed to test how effectively they recovered the electrical resistivity distribution of the subsurface. Along with implementing these schemes, we tried to quantify their accuracy, to provide robust support for the conclusions extracted from the study in this paper.

3.1. Synthetic Model

For the synthetic tests, we used some variants on an Oblique Conductor Model, which was based on the model from Ledo [46], as illustrated in Figure 1. The model comprises a regional 2-D resistivity structure formed by two half-layers of 1000 Ω·m (west side) and 100 Ω·m (east side) with their bottom at 15.75 km depth. The western half layer was covered by a horizontal layer of 500 Ω·m with its bottom at 0.5 km depth, while the eastern half layer was covered by an inclined layer of 10 Ω·m with a dip angle of approximately 3°. In the central domain, a conductive block (3 Ω·m) of approximately 20 × 4 × 4 km³ was embedded in the regional 2-D structure with its top at 3 km deep. The major axis of the conductive block obliquely crossed (45°) the contact plane of the two half layers. The bottom of the model (below the depth of 15.75 km) was set as a half space of 100 Ω·m.

The forward responses of the synthetic model were calculated for an 18 × 10 station array with a 3 km station distance that covered the model subsurface structure using the ModEM code [15,16]. The 23 periods used were evenly spaced on the logarithmic scale in a range of 0.003 and 1000 s. The central domain of the grid comprised 50 × 30 × 51 cells with a horizontal discretization of 1 km × 1 km, which was padded by five planes from the edges of the central domain outward to the boundaries with the cells lateral size increasing by a factor of two. The thickness of the top layer was 100 m, and the thicknesses of the subsequent layers increased by a factor of 1.2 in the z direction. The large 3-D grid then consisted of 60 × 40 × 51 cells in the east–west, north–south, and vertical directions, respectively.

3.2. Sensitivity Analysis

The examination scheme started with a nonlinear sensitivity test, which was performed to identify the spatial responsiveness of the oblique conductor and to characterize and quantify the effects of this conductor on the MT and MV transfer functions. The nonlinear sensitivity test involved two steps: (1) we calculated the forward responses of the background regional 2-D resistivity structure without anomaly and that of the regional 2-D resistivity structure with the embedded oblique conductor; (2) the resulting forward responses were used to calculate the differences and sensitivity values for each period and station following a formula introduced by Campanyà et al. [33].

To highlight the sensitivity of the data responses to the anomaly concerning the presence of the data errors, the differences between the forward responses were divided by the assumed error. We defined the errors for the impedance tensor data relative to the amplitude of the Z elements; namely, error levels of 0.03 × |Z_xy× Z_yx|^1/2 were adopted for all the impedance components. For the dimensionless tipper vector W responses, we adopted constant error bounds of 0.03. Although the phase tensor is a dimensionless quantity too, the error levels of the phase tensor responses here were set at 0.05 × |Φ_xx× Φ_yy|^1/2 for all Φ components instead of a constant error because the relation of the amplitude of a phase tensor element to the changes in the subsurface resistivity is highly nonlinear [37].

The contoured pseudosections of the effects of the conductive block are shown in Figure 2. Each subdiagram contains two perpendicular contoured pseudosections, the distance along the east–west profile (y = 0) on the abscissa, the distance along the south–north profile (x = 0) on the ordinate, and the period (on a decadic logarithm scale) along the vertical coordinate. The white zones denote the difference values (between −1 and 1) within the assumed errors, showing the periods for the various components at the sites that are insensitive to the anomaly. The blue and red zones represent the periods at highly sensitive sites, where the responses to the anomaly were several standard deviations outside of the assumed errors. From the calculated differences for each period and station, it can be seen that the effect of the anomalous block at a given site and period was determined by the transfer function types that were employed. Specifically, the ρ responses were similar to the Z responses, the anomaly affected a long range of periods for these two responses, and the effect continued to the longest period. For the different components of the impedance tensor Z, the period range affected by the same anomaly for on-diagonal components was narrower than that for the off-diagonal components Zo. In addition, the φ responses (impedance phase) were similar to the Φ responses; this is because the phase tensor Φ is a function of the impedance and can be written as a diagonal matrix with its components equal to the tangents of the conventional impedance phases in the layered 1-D structure [10]. Meanwhile, all of these MT responses were more sensitive at the sites above the anomaly. However, the W responses were more sensitive at the sites located at the edge of the anomaly due to the absence of a vertical field above the laterally uniform earth, where the greatest effects in the W were off to the anomaly edge [47]. Furthermore, the same anomaly affected a narrower range of periods for the Φ and W responses than that for the Z responses. The reasoning behind this is that the galvanic charges imposed on the boundaries of the conductivity contrasts generate magnetic effects and electric effects simultaneously; among them, the magnetic effects of the charges fall off rapidly with the increasing period, while the local nature of the electric fields is strong distortion from the charges on the boundaries, and the electric effects are frequency-independent and persist for longer periods in the inductive case [48]. It is worth noting that the high difference values reflected in the different data responses were mostly distributed on the west side, where there were enormous electrical differences between the embedded anomaly body and the background.

The accomplished quantification for the effect of the anomaly on the different responses illustrated in Figure 3, which comprises six subdiagrams of pseudosections, represented the sensitivity of the different transfer functions’ data, respectively. The sensitivity values of the ρ and Z responses were close, which were the highest among all the data responses, the sensitivity values of the off-diagonal impedance Zo were slightly lower than those, and the sensitivity values of the φ and W responses were approximately equal to half of the Zo responses, while the lowest sensitivity values obtained by the W responses. The increased sensitivity values of the Z responses relative to the Zo responses suggest that the use of the on-diagonal components of Z introduces the complementary information that enhances the resolution. The low sensitivity values obtained by the Φ and W responses indicate the insufficiency of the information about the absolute values of the subsurface resistivity in these dimensionless physical quantities. However, both of these are primarily sensitive to the variations of resistivity structure, which can be inferred from the distribution of those sites and periods associated with the Φ and W responses sensitive to the anomaly.

As another illustration of the sensitivity analysis, we also considered the sensitivity functions, namely the Jacobian. Six sites in an oblique line crossing the buried conductor were used to implement the sensitivity function calculations of different data responses. The used periods were set up to 1 s, 10 s, and 100 s. The results of these calculations are illustrated in Figure 4, which give the data perturbation resulting from a model resistivity perturbation. The gradients of the Zo and Z responses were the highest among all the data responses, and the gradient of the ρ + φ responses had moderate values; the gradients of the W and Φ responses, however, were three orders of magnitude lower than that of the Z responses. The calculation results suggest that the same model resistivity perturbation could result in distinct perturbations for different data responses. The low perturbations obtained by the W and Φ responses indicate the deficiency of the sensitivity to the model resistivity variations, while the medium–high perturbations obtained by the ρ + φ, Zo, and Z responses indicate sufficient sensitivity to the model resistivity variations. Although the sensitivity magnitude of the Zo was similar to that of the Z, the range of perturbations for the Z responses was a little bit larger than that for the Zo responses. This effect of the sensitivity enhancement was consistent with that of the sensitivity quantification test mentioned above (as shown in Figure 3).

3.3. Three-Dimensional Inversion Tests

For the subsequent numerical experiment, we performed a 3-D inversion of different transfer functions’ data utilizing a quasi-Newton method based on the L-BFGS formula, which was achieved by a modification to the ModEM code [16]. The solution mesh was chosen in common with the forward mesh described above, which comprised 60 × 40 × 51 cells beneath the station array. Five inversions were performed to examine and evaluate the capability of different transfer functions’ data to characterize the subsurface by themselves, that is, the inversions of the ρ + φ, Zo, Z, W, and Φ responses individually. Three additional inversions were performed to appraise the advantages and disadvantages of complementing the Zo, Z, and Φ responses with the W responses. Hence, the final inversions were (1) ρ + φ, (2) Zo, (3) Z, (4) W, (5) Φ, (6) Zo + W, (7) Z + W, and (8) Φ + W. For the subsequent inversions, all these used responses were perturbed by the same scheme described in the sensitivity analysis. The heavily penalized deviations from a prior model can lead to erroneous and misleading 3-D inversion results, particularly in environments with strong conductivity contrasts; hence, the prior model setting is important in retrieving the subsurface [28,37]. Thus, we executed a test to choose the prior model by inverting different transfer function data with four homogeneous half-space models of 10, 100, 300, and 1000 Ω·m, respectively. The results of the test are displayed in Figure S1, where the inversion results using a prior model of 100 Ω·m obtained the recovery of the electrical resistivity values most similar to the synthetic model. Based on this test, the chosen prior model for all the presented inversions was set as a homogeneous half space of 100 Ω·m. In addition, a cooling method for the damping factor was adopted to improve the stability of the inversion, and the initial damping factor value was set to 10. In addition, we inverted the Z and W responses with different model covariance smoothing parameters of 0.1, 0.2, 0.3, 0.4, 0.5, and 0.7 to discuss the influence of the model covariance matrix on the inversion results, as the model norm’s behavior is controlled by the model covariance matrix. The results of this test are displayed in Figure S2, where a model covariance smoothing parameter of 0.3 yielded the most effective inversion result. Therefore, the model covariance smoothing parameter for all the presented inversions was uniformly set to 0.3, with two smoothing passes in each direction.

The corresponding final inversion models were obtained when the acceptable overall normalized root mean square (nRMS) misfit was achieved. The fit of the data of each inversion process is shown in Figure 5; that is, the summed nRMS was calculated at each station and provided valuable information on how the data misfit was distributed across the array. It can be seen that the misfits were distributed smoothly through the array for different transfer functions’ type. Even though some larger nRMS associated with the inversion results of the ρ + φ responses were observed, the particular transfer functions’ type or station locations had no overfitting problem. These overall good fits denote that all the inversion results have a lower average variance between the original synthetic model and the final inversion model.

Table 1. Objective function value (ϕ) obtained from the inversion for different transfer functions.

Data Types	ρ + φ	Zo	Z	Φ	W	Zo + W	Z + W	Φ + W
ϕ initial	187.07	806.22	407.49	77.95	11.42	408.82	275.37	44.69
ϕ final	2.16	0.77	0.89	0.43	0.55	0.83	1.61	2.08
Iterations	76	130	93	75	31	87	72	123

shows how the objective function value changes during the iteration for each type of inversion. It can be seen that all the inversions can converge to a similar level of the objective function value, which suggests the inversions have worked properly. However, for the inversion of the W, Φ, and Φ + W responses, the decreases in the objective function values were not significant. In addition, the use of the W responses as complements to the Zo responses significantly accelerated the convergence of the inversion (i.e., reduced the maximum number of iterations from 130 to 87). However, this situation was completely opposite in the joint inversion of the two types of dimensionless quantity (i.e., the Φ + W responses).

3.3.1. Results of Inverting the MT and MV Transfer Functions Individually

The models resulting from the inversion of the MT and MV transfer functions data individually are displayed in Figure 6, where Figure 6a1–a6 illustrates the electrical resistivity model beneath a W–E profile (y = −1.5 km), and the following Figure 6b1–b6,c1–c6 represents the horizontal slices at 400 m deep and 4000 m deep, respectively. The corresponding final inversion models of all the responses characterized the main structures and resistivity contrasts of the synthetic model in approximately the correct location. The inversion results of the Zo, Z, and Φ responses recovered the closest resistivity value to the true value of the 3 Ω·m for the oblique conductive block and well constrained the lateral extent of the oblique conductive block, although the absolute resistivity of the eastern conductive block was overestimated when inverting the Zo responses. Although the inversion results of the Zo responses depicted the boundaries of the oblique conductive block in the background of high resistivity (1000 Ω·m), the boundary characterization of the oblique conductive block in the background of low resistivity (100 Ω·m) was distorted. Additionally, the inversion results of the Φ responses better defined the upper and lower boundary of the conductive block in the background of low resistivity (100 Ω·m, which was the same as that in the prior model), whereas the lower boundary of the recovered conductive block appeared at a shallower depth than in the real model when the Z responses were used in the inversion. Note that the inversion results of the Φ responses defined a wrong upper boundary in the background of high resistivity (1000 Ω·m); the reasoning behind this could be related to the initial resistivity of the prior model.

In contrast, the corresponding final inversion models of the ρ + φ and W responses could only infer the lateral extent of the oblique conductive block; although the magnitude of the resistivity contrast between the two half layers (including that between the two covers) was reproduced correctly. When only inverting the W responses, the inversion model was unable to define the continuity of the oblique conductive block; the resulting conductive block was characterized as two halves (high resistivity in the western part and low resistivity in the eastern part). This is because the W responses lacked information about the absolute resistivity values but contained information about the variations in the relative conductivity in the subsurface. There was a serious distortion in the characterization of the vertical boundary of the conductive block when inverting the ρ + φ responses, which manifested as redundant anomalies below the oblique conductive block. The reasoning behind the poorer recovery of the resistivity values and geometry of the oblique conductive block with the ρ + φ responses inversion is that the imaginary part of the Jacobi matrix was ignored when processing the apparent resistivity and impedance phase data in this paper; therefore, the accuracy of the inversion decreased. The outcome of our study indicated that the recovery of the subsurface resistivity structure showed a strong dependence on the type of transfer functions.

3.3.2. Results of Inverting the MT and MV Transfer Functions Jointly

The joint inversion results of the MT and MV transfer functions data are displayed in Figure 7 and are compared with the results of inverting the MT transfer functions data individually. It can be seen that all the joint inversion results reproduced a more accurate magnitude of the resistivity contrast between the two half layers. The use of the W responses as complements to the Zo and Z responses in the inversion slightly improved the resolution of the oblique conductive block. In these cases, the recovery of the resistivity values and geometry of the oblique conductive block was significantly improved. Particularly for complementing the Zo responses with the W responses, a more homogeneous resistivity distribution of the oblique conductive block was defined. However, the lower boundary of the oblique conductive block seemed to appear at a shallower depth when complementing the Φ responses with the W responses. These results indicate that the W responses were insufficient to characterize the variations of conductivity in the vertical direction, but the corresponding final inversion models preferably constrained the lateral variations of the conductivity when inverting the W responses as a complement to the Zo, Z, and Φ responses.

3.3.3. Assessment of the Model’s Accuracy

In order to quantify how well the models obtained from the inversion process reproduced the synthetic model, we calculated the differences between them in model space. The results with a focus on the vertical structures below the W–E profile and the lateral structures at a depth of 4000 m are shown in Figure 8. The differences between the synthetic model and the inversion models were calculated using the decadic logarithm electrical resistivity values. The colors indicate the resistivity biases between the inversion model and the synthetic model, where the red and blue colors, respectively, indicate the results that were more conductive and more resistive than the synthetic model, respectively, and the white color indicates an almost-perfect fit between the final inversion model and the original synthetic model. The average differences for each type of inversion are shown in

Table 2. The average differences between the final inversion models and the original synthetic model.

Data Types	Zo	Z	Φ	Zo + W	Z + W	Φ + W
Difference	0.127	0.139	0.124	0.119	0.132	0.172

. It can be seen that all the inversion models better defined the geometry of the conductive block, although the electrical resistivity values of the conductive block edges were more resistive than in the synthetic model. The results from the inversion processes complementing the W responses with the Zo, Z, and Φ responses significantly improved the recovery of the electrical resistivity values of the regional 2-D resistivity background, showing a more homogeneous resistivity distribution and correctly reproducing the magnitude of the resistivity contrast between the two half layers. However, the electrical resistivity values associated with the lower boundary of the conductive block were more resistive than in the synthetic model when the W responses were involved in the inversion of the Φ responses.

4. Discussion

In this work, we performed sensitivity tests on the different MT and MV transfer functions, including the ρ, φ, Zo, Z, W, and Φ responses. The results of our numerical experiment suggest that the difference in the sensitivity of the Z and Zo responses to the anomalies directly leads to the difference in their capabilities to recover the subsurface electrical resistivity. The capacity for recovering the detail of the geoelectrical structures (Figure 6) was remarkably improved when the on-diagonal components of Z were used. Although the sensitivity of the ρ responses was similar to that of the Z responses, the ρ + φ responses were not good types of data to define the boundaries of the subsurface geoelectrical structure and even resulted in redundant artifacts under the oblique conductive block (Figure 6), which may be caused by ignoring the imaginary part of the Jacobi matrix. Inspecting the sensitivities obtained from the W responses, the sensitivity associated with the oblique conductive block was the lowest, and the locations and periods responsive to the anomaly body appeared as a narrow zone. This property of the W responses hence resulted in the worst recovery ability of the subsurface geoelectrical structure. For the Φ responses, the sensitivity in connection with the anomaly body was moderate, and the locations and periods responsive to the anomaly body were sensitive at the edges of the structures, making the Φ responses a good type of data to define the boundaries of the structures. The consequences of this property in the inversion process can be seen in Figure 6, Figure 7 and Figure 8, where the models obtained by inverting the Φ responses were closer to the original synthetic model; in particular, the boundary of the oblique conductive block was better than that by inverting other data types. Although the W responses were intrinsically less sensitive to the existence of the anomaly than the Zo, Z, and Φ responses, it effectively improved the ability to recover the detail of the geoelectrical structures when it was complemented with the Zo or Z responses. This improvement was more remarkable for the off-diagonal components of Z, the models including the W responses obtained the horizontal boundary of the oblique conductive block closer to the original model (Figure 7), and the difference between the final inversion model and the original synthetic model was smaller than that obtained by inverting the Zo responses individually (Figure 8 and Table 2). However, the use of the W responses as an accompaniment to the commonly used Φ responses obtained the incorrect bottom boundaries of the geoelectrical structures, which increased the difference between the final inversion model and the original synthetic model (Figure 8 and Table 2), even though the recovery of the background resistivity was better.

Many previous studies have emphasized that the use of diagonal elements in the inversion of the full impedance tensor Z possesses a superior capacity to reproduce the dimensionality of the 3-D geoelectrical structures [28,29]. These conclusions can be validated by the results from the inversion of the Z and Zo responses in our work; it can be seen that the use of the Z responses aided the superior determination of the conductive block than when using the Zo responses, including the characterization of the geometry and the absolute electrical resistivity values (Figure 6). The phase tensor Φ was more sensitive to the vertical contrasts and was unaffected by galvanic distortion; the approach of inverting the phase tensor Φ is highly recommended for 3-D inversion [37,39]. Although the results of inverting only the Φ responses presented a stronger dependence on the resistivity level of the prior model (Figure S1), complementing the Φ responses with the W responses can lead to the superior recovery of the absolute resistivity structure [38]; this inference is supported by the recovery of the background resistivity shown in our work (Figure 7 and Figure 8). The tipper vector W was more indicative of the lateral contrasts; thus, the results of inverting the W responses as a complement to the Zo, Z, and Φ responses provided preferable lateral contrasts of conductivity, particularly for the recovery of the background resistivity model (Figure 7 and Figure 8).

5. Conclusions

This study provided an appraisal of the MT and MV transfer functions’ selection in a 3-D inversion by utilizing a quasi-Newton method based on the L-BFGS formula. The inversion result showed that the 3-D inversion of either the Z or Φ responses presented a superior ability to recover the electrical resistivity values of the oblique conductive block as well as the background resistivity model. Moreover, the joint inversions of the Zo, Z, or Φ responses with the W responses possessed a superior imaging of the geometry and horizontal continuity of the oblique conductive block. From these results, we conclude that the ability to recover the subsurface electrical resistivity depends on the sensitivity of the transfer function. The greater the transfer function resistivity perturbation produced by the model resistivity perturbation, the stronger the ability of the transfer function to recover the underground structure. Comparing the effect of the individual and joint inversion of different transfer functions is an important way to determine the advantages of each transfer function and achieve a more reasonable interpretation.