Parallel Simulation of Audio- and Radio-Magnetotelluric Data

: This paper presents a novel numerical method for simulation controlled-source audio-magnetotellurics (CSAMT) and radio-magnetotellurics (CSRMT) data. These methods are widely used in mineral exploration. Interpretation of the CSAMT and CSRMT data collected over an area with the complex geology requires application of effective methods of numerical modeling capable to represent the geoelectrical model of a deposit well. In this paper, we considered an approach to 3D electromagnetic (EM) modeling based on new types of preconditioned iterative solvers for finite-difference (FD) EM simulation. The first preconditioner used fast direct inversion of the layered Earth FD matrix (Green’s function preconditioner). The other combined the first with a contraction operator transformation. To illustrate the effectiveness of the developed numerical modeling methods, a 3D resistivity model of Aleksandrovka study area in Kaluga Region, Russia, was prepared based on drilling data, AMT, and a detailed CSRMT survey. We conducted parallel EM simulation of the full CSRMT survey. Our results indicated that the developed methods can be effectively used for modeling EM responses over a realistic complex geoelectrical model for a controlled source EM survey with hundreds of receiver stations. The contraction-operator preconditioner outperformed the Green’s function preconditioner by factor of 7–10, both with respect to run-time and iteration count, and even more at higher frequencies.


Introduction
Controlled-source audio-magnetotellurics (CSAMT) is a popular method for mineral exploration [1][2][3][4][5]. In principle, a non-uniform plane wave at a study area is generated by two orthogonal transmitter lines located at a sufficient distance (in the far-field zone). A typical operation frequency range is 0.1-10 kHz, which is roughly the high end of the quasi-stationary regime. Transfer functions between electric and magnetic field components, Ex, Ey, Hx, Hy, and Hz are computed from measured data and then inverted in order to reconstruct the subsurface electric conductivity. CSAMT is known to provide better data quality in the AMT dead bands (0.1-5 Hz and 700-3000 Hz) than conventional AMT surveys and it requires relatively inexpensive transmitter operations [2,5,6].
Closely connected to CSAMT is the controlled-source radio-magnetotellurics (CSRMT). This method is gaining popularity in the near-surface geophysics [7][8][9][10][11]. It operates at much higher frequencies (typically, up to 250 kHz-1 MHz) and thus provides excellent spatial resolution. The biggest difference with CSAMT is the displacement currents in the air [10,12] cannot be neglected to In CSAMT/CSRMT, the five components of the electromagnetic filed are measured (Ex, Ey, Hx, Hy, Hz) and, after preprocessing, converted to the components of the impedance and tipper. The off-diagonal components of the impedance are employed most commonly, that is [2], = , where subscripts 1 and 2 refer to one of the two orthogonal transmitter lines. Thus, computation of the components of the impedance or tipper requires simulation of the electric and magnetic fields and at several locations and frequencies. We consider electromagnetic modeling in 3D heterogeneous isotropic media in the frequency domain. Assuming time dependence of , the electric field, ( , , ), satisfies the following system of partial differential equations, where is the source angular frequency, is the magnetic permeability of the free space, ( , , ) is the conductivity model, and ( , , ) is the source current density. This is a linear system of three scalar equations with respect to three scalar entries of ( , , ). We solve this partial-differential equation system in a bounded domain, completed by zero Dirichlet boundary conditions. The magnetic field is then computed via Faraday's law, Since geological models predominately vary in the vertical direction, we can always split the conductivity model into background and anomalous parts, We further assume that the following double inequality holds, which controls the contrast of anomalous conductivity with respect to the background. Given a non-uniform computational grid with cells, we apply the conventional edge-based second-order FD discretization. The FD method approximates the unknown electric field = ( , , ) with a finite set of discrete values { , , } [26]. Each discrete value is attached to the respective edge of the FD grid, Figure 1. The actual discrete equations can be found for example in [27]. Let us form the unknown vector of the discrete electric fields { , , } introduced above. Now we can write the FD discretization of (2) in a matrix form: where is a complex symmetric system matrix with at most 13 nonzero entries per row, corresponding to the total conductivity distribution. We denote the size of this system as , ≈ 3 . Let be the FD system matrix, corresponding to the background conductivity model. Importantly, this matrix can be implicitly factorized using a fast direct algorithm, and the action of the inverse matrix can be efficiently computed [23,28]. As a result, it can be used as a preconditioner to (6): We will refer as the Green's function (GF) preconditioner. The complexity of applying the GF preconditioner is ⁄ , and auxiliary memory required is near 3 only. Applying the analysis presented in [23] the condition number of the GF preconditioned system can be estimated as follows, This result implies that convergence of an iterative solver applied to (7) has minor or no dependence on the grid size and cell aspect ratio, as well as the frequency, while it degrades on models with high-contrast bodies.
To minimize the impact of high-contrast bodies, another preconditioner could be constructed. Denote as and diagonal matrices, corresponding to full and background conductivity, respectively. Let us define the modified FD Green's operator and diagonal matrices according to the following formulae: Using this operator, Equation (7) can be written in an equivalent form as follows (see Appendix A): By introducing a new operator, = , we rewrite Equation (11) as follows: We will refer this system a contraction operator (CO) preconditioned system, since it was shown in [23] that ‖ ‖ < 1. Moreover, it can be proved that this preconditioned system has a smaller or equal spectral condition number than that of the GF preconditioner, implying faster or equal convergence of the iterative solvers, The complexity of applying the CO preconditioner is ⁄ as well. The preconditioners in Equations (7) and (12) were incorporated into the BiCGStab iterative solver [29]. We used the complex-valued version of the solver with the standard complex-valued dot product.
Practical controlled-source modeling requires excessive gridding near the source location. To avoid this, we preferred secondary field modeling, i.e., the electric field ( , , ) is represented as a sum ( , , ) + ( , , ), where ( , , ) is the response due to background conductivity model known analytically [30,31]. In this case, the secondary field ( , , ) will be the unknown function and its FD discretization is performed. The actual source ( , , ) in Equation (2) is substituted with the secondary source ( , , ) ( , , ),

Applicability of Quasi-Statinary Simulation
As the first step, we tested applicability of the quasi-stationary simulation in the context of controlled-source electromagnetics. The displacement current in the air must be considered when the electromagnetic field generated by a high-frequency dipole oscillator is measured at large offsets [10,12]. Otherwise observed components of the electromagnetic field cannot be matched to the computed quasi-stationary ones; this is known as the propagation effect. Formally, it is achieved by replacing term in Equation (1) with − , where is the vacuum dielectric permittivity and ≥ 1 is the relative dielectric permittivity. It complicates the 3D numerical simulation considerably. However, in contrast to the individual components, the surface displacement-current impedance was almost identical to the quasi-stationary one, at least in a 1D Earth. To illustrate this point, we considered an X-oriented horizontal electric dipole (HED) located at the origin of Cartesian coordinate system and two receiver stations ( Figure 2). Figure 2. Impact of the displacement currents on the electromagnetic field and the surface impedance on the top of a layered Earth. The left column of panels illustrates computations relative to a short transmitter-receiver separation (receiver Rx-1). The right column of panels presents computations for a long transmitter-receiver separation (receiver Rx-2). Panels (a,b) depict geometry of the simulation, the top view. An X-oriented horizontal electric dipole (HED) was placed at the origin. Receiver station Rx-1 was located at (450, 750) m, whereas receiver station Rx-2 had coordinates (450, 4500) m. Panels (c,d) depict magnitudes of individual components of the electric field. Panels (e,f) depict Cagniard's apparent resistivity ( ). Panels (g,h) show impedance phase. In the legends, QS stands for the quasi-stationary computations.
The first station corresponded to relatively short offset (X = 450 m, Y = 750 m) and the second one imitated a typical CSAMT offset (X = 450 m, Y = 4500 m). We computed the electromagnetic field for the two cases. In the first case, the model was the one in Table 1. Computations were conducted in the quasi-stationary regime, where the squared wavenumber of -th layer was given by . The air conductivity was set to =10 S/m. In the second case, the squared wavenumbers of each -th layer were defined as − . We set = 1 in the air, and = 4, = 1,2, … in the Earth. The computations were performed by the 1D code of [13]. Computed curves are presented in Figure 2. At frequencies higher than 7 kHz (the close receiver) and 27 kHz (the distant receiver), the quasi-stationary electromagnetic field (dashed lines) differed considerably from the displacement-current field (solid lines). However, ratios of the horizontal component remained essentially the same.

Implementation
Forward modeling of a single source involves the following steps: • Modeling grid preparation, • Resistivity model resampling, • Right-hand side computation, • Secondary electric field computation, and • Computation of the total electric field, recovery of the magnetic field.
Our implementation was in C++. The most computationally demanding steps are the right-hand side (RHS) and secondary field computation. To reduce the computational burden, these steps were parallelized using OpenMP shared memory implementation. Simulation of multiple sources and frequencies was parallelized using the Message Passing Interface MPI. In [32], we demonstrated good scalability of this scheme.

Code Verification on Simple Models
Our code in general has been validated and benchmarked against analytical layered-earth plane-wave solutions in [23]. Verification of the code in the low-frequency controlled-source electromagnetics scenario was performed in [32]. Here, we conducted several tests concerning particularities of the CSAMT/CSRMT, mainly high frequencies and high conductivity contrast.
We used Consortium for Electromagnetic Modeling and Inversion's software PIE3D for benchmarking. PIE3D was rigorously verified and has been used in production for years [25,33,34].
In this test, we used a resistivity model consisting of a 1D background model and a 100 Ωm parallelepiped body. The background model is presented in Table 1   The transmitter was an X-orientated point electric dipole located at (32.54, −553.5) m. The receiver was at (200, 80) m. In this test, we used seven frequencies: 0.1, 0.2, 0.5, 1, 2, 5, and 10 kHz.
Finite-difference computations were performed with large and fine numerical grids. The core domain was the same at each frequency 478 × 158 × 150 m 3 . It included both the anomalous body and receiver location (Figure 4). PIE3D ran with grid step size h = 2, meaning that the body was discretized into 66 × 28 × 16 cubical cells. Computed electromagnetic responses are compared in Figure 5. . Benchmark against PIE3D software. "IE" and "FD" mean the integral-equation software (PIE3D) and the finite-difference software (our code), respectively. "Re" and "Im" mean the real and imaginary parts, respectively. Note that the difference curves are magnified by factor 10.
We observed a close match between individual components. The normalized misfit between IE and FD responses, that is ‖ − ‖/‖ ‖, varied from 0.2% to 0.6%. This is encouraging, taking in mind that the computations were done by two very different approaches with no shared code.

Conductivity Model of Aleksadrovka
We evaluated our numerical code on a realistic 3D model using an acquisition geometry from a real geophysical survey. The model mimicked composition of Moscow State University's geophysical test camp near the village of Aleksandrovka in Kaluga Region, Russia. The acquisition geometry is presented in (Figure 6). In this work, we did not compare modelled versus measured data. Figure 6. Map of the study area from which the acquisition geometry was taken for the numerical simulation. The boxes depict receiver stations. The receiver lines are numbered from 1 to 8.
The receiver grid consisted of 192 receiver stations and two orthogonal transmitter lines (see Figures 6 and 7). All receivers had azimuth 68.5°. Coordinates of receiver stations and transmitter lines can be found in the supplementary materials File S1. The five components of the electromagnetic field were simulated. The model used in this study was compiled from various geologic studies, geophysical surveys, including seismics, and drilling data. Aleksandrovka area is located in Moscow syneclise and has relatively simple morphology of pre-Quaternary layers, but the shallow part of the subsurface has a more complicated structure due to Moscow glaciation [35]. The reference 1D geoelectrical model is based on a 300 m borehole drilled in the camp for water supply purposes. The lithology column is shown in Figure 8. The final 1D reference model is presented in Table 2. The 1D horizontally layered reference model was deformed based on the topography of the top of the first layer of resistive carbonates. This reference surface was clearly identified from previous CSRMT survey. The most prominent feature of it is the trough oriented along the NW-SE direction, which has the maximal relative depression of about 60 m ( Figure 8). Other layers with thicknesses taken from the 1D model conform to the reference surface.
Finally, we added to the model inhomogeneities that represent composition of the shallow part of the subsurface. The subsurface has a relatively complicated geoelectrical structure because of Moscow glaciation. This part has been characterized previously with a detailed CSRMT survey in the frequency range 5-1000 kHz (the far-field responses only) followed by a conventional 2D plane-wave inversion. Several superficial high-resistive sand bodies were delineated on the interpreted 2D cross-sections ( Figure 9).
The biggest body was composed of Quaternary glacial or alluvial sands of 160 Ωm with its top at a depth of 20 m. The eastern part of it has an elongated shape, with approximately 100 m along its short axis and thickness of 20 m. This body was inserted into the 3D model. Finally, a 6-m-thick layer of 19 Ωm was introduced across the top of the 3D resistivity model. The final 3D model is given in Figure 10a,b together with receivers and transmitters. We believe it represents all main features of the Aleksandrovka site. The conductivity model in the Visualization Toolkit VTK legacy file format (https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf) is available in the supplementary materials File S1 attached to this article.

Numerical Simulation of Aleksadrovka
We compared efficiency of the sequential 3D modelling code both for GF (7) and CO (12) preconditioners. The transmitter was a point dipole located at the southern end of transmitter line Tx-2 (see Figure 3), with coordinates (1091.101, 812.101, 0.002) and azimuth 155.6°. We selected three frequencies, 192, 320, and 576 Hz. Parameters of the FD grids are presented in Table 3. The code was running on a single thread of a compute node equipped with double 12-core Intel Xeon E5-2680v3 processors running at 2.5 GHz and 128 GB of memory. The relative accuracy of the BiCGStab was set to 10 −10 . Performance of the two preconditioned iterative solvers is presented in Table 4. The table also includes CPU time needed to prepare the RHS. We observed an acceleration of 7×, 9×, and 10× at the three selected frequencies, respectively, both with respect to the iteration count and wall-clock time. Computational load per iteration for the CO preconditioner was only 7-10% larger compared to the GF one. Thus, a reduced number of iterations almost directly translated to the shorter compute time. The GF solver did not converge at 576 Hz frequency for 5000 iterations, which was set as the limit. We attributed slow converge of the GF solver to the highly conductive 1.5 Ωm layer which produced a strong horizontal resistivity contrast in the model. It should be mentioned that topography slowed the computations considerably since performance of the GF preconditioner deteriorated most severely when the airground interface was not a horizontal plane. In this case, the time gain of the CO over GF preconditioner was expected to be even greater than one observed in this test.
It is interesting that the time of computing the right-hand side dominated by computation of the background electric field was comparable to the time spent on the solution of a system of linear equations. This is one of the disadvantages of the secondary-field approach. However, RHS computation was easily reduced by using shorted Hankel transform filters and multithreading (see Section 2.2). In the next tests, we applied the CO preconditioned iterative solver only.
In the second set of computations, the whole survey was simulated. The electromagnetic field was computed at 192 receivers. We used 18 frequencies with range from 192 Hz to 550 kHz. It should be emphasized that the electromagnetic field above 25 kHz was not quasi-stationary for our setup, but impedance computed from individual components still can be used to interpret the real measurements (see [13]). Two transmitter lines were used to compute the two polarizations. During computations, the lines were broken into 26 straight segments (15 for Tx-1 and 11 for Tx-2), which had lengths between 40 and 185 m, and then numerically integrated along the transmitter lines. Thus, the total number of individual forward problems was 468. Parameters of the FD grids at different frequencies are given in Table 5. Simulation was conducted on 26 nodes equipped with double 12-core Intel Xeon E5-2680v3 (24 cores per node) processors running at 2.5 GHz, 128 GB RAM, and InfiniBand interconnect. Each compute node was processing the 18 forward problems in the frequency range 192 Hz-550 kHz. There was message passing between nodes during computations. At each node, all 24 cores were working on each forward problem using OpenMP parallelization. The overall computations lasted for 24.5 h. Duration of the individual computation phases is given in Table 6. We observed that computation of the right-hand side (essentially, the background electric field) was comparable to the time allocated to the solution of a linear system. At frequencies above 150 kHz, the right-hand side dominated the other computations. It was caused by a combination of large FD grids and excellent performance of the CO preconditioned solver at the higher frequencies, at which core domains were very thin and enclosed within the top part of the model with moderate conductivity contrast. It is, however, of minor concern for inverse problems, because the background electromagnetic filed can be stored on disk. The cost of its computation is amortized across many forward problems solved over the course of the inverse problem.
Apparent resistivity and impedance phase along receive line No. 6 are presented in Figure 11. The two polarizations resulted in very similar cross-sections. Apparent resistivity approached the limit of 18 Ωm at frequencies above 25 kHz, indicating the high-frequency electromagnetic field did not penetrate below the top layer. Figure 11. Results of the numerical simulation. Pseudo cross-sections along receiver line No. 6. Top two panels: XY and YX apparent resistivity. Bottom two panels: Impedance phase components, XY and YX.
Maps of the XY apparent resistivity and impedance phase are presented in Figure 12. The map at 192 Hz had a clear correlation with the topography of the first limestone layer (Figure 8), whereas the 1.5 Hz data clearly indicated the high-resistivity sand body (Figure 10b).

Discussion
In this paper, we considered an effective approach to modeling the data acquired by controlled-source audio-magnetotellurics (CSAMT) and radio-magnetotellurics (CSRMT) methods, which are widely used in mineral exploration. We assessed performance of two preconditioned iterative solvers for CSAMT/CSRMT simulation introduced earlier in [21]. Further, we presented results of numerical modeling of a real CSRMT survey performed near Aleksandrovka, Kaluga Region, Russia. The 3D model we used was based on 2D and 1D inversion as well as drilling data.
Collectively, this study demonstrated that the CO preconditioner is highly efficient even in extreme conditions of high excitation frequencies and large conductivity contrasts, encountered in CSAMT/CSRMT surveys. In contrast, the GF preconditioner degraded severely at the higher frequencies. This observation has not been published earlier as far as the authors are concerned.
From our experience, extremely large FD grids (above 50 M discrete unknowns) are easily encountered at high frequencies typical for CSRMT (above 10-50 kHz, depending on the subsurface resistivity). However, these conditions were pathological in the sense that the quasi-stationary regime is not valid anymore, and the displacement currents must be considered. The only reason such computations make sense is the quasi-stationary impedance was similar to the displacement-current one, at least in the 1D Earth, as demonstrated by [13].
From the geophysical standpoint, the data at frequencies above, perhaps, 25 kHz were redundant, because the 3D conductivity model did not contain very shallow inhomogeneities (<10 m). Our purpose was to evaluate numerical grids and simulation time if such high frequencies will be needed in the future.
Finally, we emphasize that the results of our study apply to other geophysical methods that require quasi-stationary electromagnetic modelling, specifically MT, AMT, CSMT, CSAMT, CSEM, borehole methods, and others.
Funding: This research was funded by Russian Science Foundation, grant number No. 16-11-10188. Acknowledgments: The authors thank M.S. Zhdanov for the opportunity to benchmark our code with PIE3D software, and M. Čuma for their help in performing such calculations. This work has been carried out using computing resources of the federal collective usage center Complex for Simulation and Data Processing for Mega-science Facilities at NRC "Kurchatov Institute", http://ckp.nrcki.ru/. We would like to acknowledge the Skoltech CDISE's high-performance computing cluster, Zhores [36], for providing the computing resources that have contributed to the results reported herein.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In this appendix, we show equivalence of Equations (7) and (12). Before moving on, denote = − which is a discrete secondary conductivity. We will start with Equation (12). After switching from to and multiplying by from the left, we receive, Next, we substitute and and employ commutativity of diagonal matrices, Now, we substitute and rewrite the matrices, It is easy to note that the term cancels out, and we receive, After factoring out in the left-hand side and noting that = − , we obtain which is exactly (7).