New 3-D Combined Inversion Scheme Using Response Functions Free From Galvanic Distortion

The combined inversion using distortion-free response functions is an effective 10 approach to robustly estimate the 3-D electrical resistivity structure against the 11 distortions caused by near-surface resistivity anomalies. However, previous combined 12 inversion analyses have presented a significant dependency of the inversion results on 13 initial and prior models. Therefore, in this study, we evaluated the effectiveness of the following two new types of 3-D combined inversion using distortion-free response 15 functions: one uses the phase tensor and the vertical and inter-station horizontal 16 magnetic transfer functions, while the other uses the Network-MT response functions, 17 in addition to the former. Because long dipoles are used, the Network-MT response 18 function is negligibly affected by galvanic distortion. To access the combined inversion 19 approach, we developed a novel 3-D inversion scheme combining the response 20 functions of the usual magnetotelluric measurements and the Network-MT response 21 function. The synthetic inversion analysis demonstrated that both of the proposed combined inversions can recover the characteristic resistivity distributions of the target model without a significant dependence on the initial models, at least in the shallow part. These results demonstrate that the combined inversions using only distortion-free 25 response functions have the potential to estimate subsurface resistivity more robustly 26 than what was previously thought. Furthermore, we confirmed that the combined 27 inversion using the Network-MT response function can make the resultant resistivity 28 structure closer to the actual one and enhance the stability of the inversion. This result 29 suggests that the combined use of the Network-MT response function is the preferred 30 approach.

part. These results demonstrate that the combined inversions using only distortion-free 25 response functions have the potential to estimate subsurface resistivity more robustly 26 than what was previously thought. Furthermore, we confirmed that the combined 27 inversion using the Network-MT response function can make the resultant resistivity 28 structure closer to the actual one and enhance the stability of the inversion. This result 29

Introduction 37
The three-dimensional inversion of magnetotelluric (MT) data has become common in 38 recent decades. This is because of the increase in computational resources and the 39 widespread use of practical 3-D MT inversion codes, such as WSINV3DMT (Siripunvaraporn et al. 2005; Siripunvaraporn and Egbert 2009) and ModEM (Egbert 41 and Kelbert 2012; Kelbert et al. 2014). However, the galvanic distortion of the MT 42 impedance tensor has long been recognized as a major problem in estimating subsurface 43 electrical resistivity structures. Galvanic distortion, which is caused by lateral small-44 scale inhomogeneities near the Earth's surface, distorts the electrical field and, hence, 45 the estimated impedance tensor (e.g., Ogawa 2002; Simpson and Bahr 2005). The previous studies using distortion-free response functions presented a deficiency in 65 robustness. Patro et al. (2013) and Tietze et al. (2015) showed that the resultant 66 electrical resistivity structure of the inversion using either or both PT and VMTF 67 depends significantly on the initial and prior models. They concluded that a sensible 68 selection of the prior and initial models is essential for recovering a reliable resistivity 69 structure. Campanyà et al. (2016) demonstrated that the inversion only using VMTF and 70 HMTF failed to recover the subsurface resistivity variation in the vertical direction and 71 the absolute resistivity values. They suggested that the inversion using either VMTF or HMTF is reliable only when additional information defining the background resistivity 73 structure is provided. As discussed in the previous studies, the dependency on the prior 74 and initial models should be caused by the low sensitivity of PT, VMTF, and HMTF to 75 the absolute values of the subsurface electrical resistivity. In addition, the combined use of NMTRF is advantageous because it enhances the 94 sensitivity to deep structures. Long dipole measurements of the NMT method make the 95 signal-to-noise ratio of the observed electrical field higher than that of the standard MT 96 measurement. In addition, well-maintained telephone lines enable us to perform long-97 term observations from several months to several years 98 Uyeshima et al. 1989), allowing us to increase the signal-to-noise ratio by stacking the 99 data. Consequently, we can estimate NMTRF with small errors, including periods 100 longer than several thousand seconds, which assists in estimating a reliable electrical In the following sections, we first present the algorithm of the combined inversion 116 scheme used in this study, including the treatments of NMTRF in the inversion. Next, 117 we describe the synthetic inversion tests to evaluate the performance of the combined 118 inversion approaches and discuss the results of the combined inversions. 119

Algorithm of the combined inversion scheme 121
We used the inversion scheme based on the 3-D MT inversion code of Usui (2015) and 122 where ∈ ℝ 3 and E ∈ ℂ denote the vector basis function (Jin 2002) and the 142 tangential component of the electric field on the -th edge of the element, respectively. 143 As shown in Usui (2015), we can obtain a linear equation for each element by using the 144 Galerkin method as follows: 145 where = (E 1 , … , E 6 ) ∈ ℂ 6 is the solution vector, ∈ ℂ 6 denotes the right-146 hand-side vector including the boundary condition, and the -th row and -th column of 147 the coefficient matrix ∈ ℂ 6×6 is expressed as follows: 148 In Equation (4), V denotes the volume occupied by the -th element. By assembling 149 all the element equations, we can obtain the following linear equation. 150 The dimension of (5) is the number of unknown tangential components of the electric 151 field at the element edges, and the coefficient matrix is a symmetrical sparse 152 complex matrix. MKL PARDISO is used to solve the linear equation. To calculate the 153 response functions as indicated below, forward computation is performed twice for each 154 frequency by appending a source field in the x-direction (E x -polarization) and y-155 direction (E y -polarization). After obtaining the electric field by solving the linear 156 equation above, the magnetic field is calculated from the electric field using 157 Faraday's law. 158 From the results of the E x -polarization and E y -polarization, we can calculate the 159 frequency response function at each observation station. VMTF ∈ ℂ 2 is calculated 160 as follows: Meanwhile, HMTF ℎ ∈ ℂ 2×2 is calculated from the horizontal magnetic field 162 components at an observation station and those at a reference station as follows: 163 where the characters with subscript r denote the magnetic field components at the 164 reference station. Further, PT ∈ ℝ 2×2 is calculated from the real and imaginary parts 165 of the impedance tensor as follows: 166 where Zxx, Zxy, Zyx, and Zyy are the components of the impedance tensor ∈ ℂ 2×2 , which 167 is computed as follows: 168 NMTRF = (Y x , Y y ) ∈ ℂ 2 is defined as the frequency response function between the 169 voltage difference V ∈ ℂ for the dipole of an NMT station and the horizontal magnetic 170 field components at a reference station (Uyeshima et al. 2001).
In the developed combined inversion scheme, voltage differences were computed using 172 an algorithm similar to that proposed by Siripunvaraporn et al. (2004). Because long 173 dipoles of the NMT method generally extend over multiple finite elements of the 174 computational mesh, dipoles are divided into segments by element edges (Figure 1). 175 Each segment is located within an element surface, which is triangular. The voltage 176 difference along a dipole can be obtained by summing the voltage differences of the 177 small segments. According to Faraday's law with the e -iωt time dependence, the voltage 178 difference dV ∈ ℂ along each segment (illustrated as a red dashed line in Figure 1) can 179 be calculated as follows: 180 where E 1 and E 2 are the electric field components along the triangle edges, L 1 and L 2 181 are the lengths of the edges, i is the imaginary unit, H n is the magnetic field 182 component normal to the triangle, and S is the area of the triangle. The magnetic field 183 component normal to triangle H n is calculated from the electric field components of the 184 element edges on the Earth's surface using Faraday's law (Equation (6)). Based on the voltage differences and the horizontal magnetic field of the two polarizations, NMTRF 186 can be calculated as follows: 187 In the combined inversion, we find the model parameters (subsurface electrical 188 resistivities) that minimize the objective function ϕ( ) ∈ ℝ: 189 where is the diagonal matrix of the diagonals that are reciprocals of the errors in the 190 where M is the number of model parameters; (17) 218

Synthetic inversion analysis and discussion 219
We performed inversions using two different datasets. The first dataset contains PT, 220 VMTF, and HMTF, whereas the second dataset contains PT, VMTF, HMTF, and 221 NMTRF. Hereon, we refer to the former and latter cases as Comb-A and Comb-B, 222 respectively. We performed a synthetic inversion analysis to assess whether the two new 223 types of the combined inversion approaches can reproduce the appropriate resistivity Meanwhile, we calculated NMTRF for ten logarithmically spaced periods from 10 to 238 10,000 s because it is usually difficult to obtain the response functions at short periods 239 at NMT stations. Because NMTRF can be estimated with small errors at periods longer 240 than several thousand seconds, as indicated in the Introduction section, we set the 241 maximum period for NMTRF to be longer than the other response functions. Because it 242 is difficult to incorporate a thin 0.1 km near-surface layer (100 Ωm) of the OC model 243 using a tetrahedral mesh, which requires a huge number of elements to represent such a thin layer, we used a forward program using hexahedral mesh (Usui 2020; Usui et al. 245 2020) to calculate the synthetic response functions. The algorithm for the forward 246 calculation using the hexahedral mesh is the same as that described in Section 2. By 247 using different meshes, we can also avoid the so-called "inverse crime." Before 248 calculating PT from the impedance tensor, we provided a galvanic distortion to the The standard deviations for NMTRF were 3% of their absolute values; that is, 0.03|Y x | 259 or 0.03|Y y |. For PT, we added Gaussian noise to the distorted impedance tensor and then calculated the PT errors by error propagation, as indicated in Patro et al. (2013). 261 The standard deviation of the noise to the impedance tensor was 3% for each 262 component, along with a floor of 3% of |Z xy Z yx | 1/2 for diagonals, following Tietze et 263 al. (2015). We started the combined inversion from three different uniform half-space 264 models (10, 100, and 1,000 Ωm) to investigate the dependency of the inversion results 265 on the initial model. In all cases, the tradeoff parameter α was set to 10.  aforementioned results suggest that the combined inversion using only the response 386 function free of galvanic distortion has the potential to robustly recover the target 387 resistivity structure, especially by using NMTRF. Therefore, we expect that the developed combined inversion approach will help in robustly recovering resistivity 389 structures by avoiding the distortions caused by near-surface anomalies.

Availability of data and materials 402
The datasets used and/or analyzed during the current study are available 403 from the corresponding author on reasonable request.      Electrical resistivity variation with a depth at (a) (x, y) = (10 km, -10 km) and (b) (x, y) = (-10 km, 10 km).
The thick black lines indicate the resistivity variation of the true resistivity structure. The broken colored and solid-colored lines indicate the resistivity variations of the resistivity structure obtained by Comb-A and Comb-B, respectively.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. supportinginformation.docx