Parallel Simulation of Audio- and Radio-Magnetotelluric Data

Abstract

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.

1. 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, E_x, E_y, H_x, H_y, and H_z 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 accurately compute the electric and magnetic fields, which complicates the numerical simulation dramatically. The authors of [13] demonstrated, however, that the full impedance and quasi-stationary impedances are almost identical in the 1D Earth. It suggests that the CSRMT measurements can be interpreted using an inversion software with the quasi-stationary forward simulation engine (note, though, [14,15]).

Numerical forward modeling plays a paramount role in both electromagnetic (EM) survey design and data interpretation. This demand has been driving the investigation of efficient modeling methods for nearly 60 years. Today, 1D and 2D frequency-domain modeling is relatively established, while research on 3D modeling algorithms is a hot topic.

The common problem in 3D modeling of CSAMT and CSRMT data is the solution of a large system of linear equations. This problem is present in all the discretization methods: Finite difference (FD), finite element (FE), integral equation (IE), and others; although in FD and FE modeling the system matrix is sparse, while in IE modeling the system matrix is dense. Consequently, design and implementation of an economical, scalable, and robust solver for this system is the key to successful 3D modeling and inversion.

At the frequencies from 0 to 1 MHz, the EM fields are mainly diffusive. Thus, the errors introduced by resistivity model resampling are local. Consequently, adaptation of the grid to a particular interface is a minor concern. We thus consider the second-order edge-based finite-difference discretization.

Computer simulation of both CSAMT and CSRMT data poses great computational challenges. The land surveys are frequently characterized by an exceptionally large resistivity contrast of the rocks, e.g., 1000 Ωm for granite and 0.01 Ωm for nickel ore [16], leading to very fine numerical grids needed to resolve the shortest skin depth. At the same time, the computational domain should be extended to include attenuation dictated by the longest skin depth, especially at frequencies above ~1 kHz, at which considerable amount of the electromagnetic field travels in the air and resistive rocks with little attenuation. Furthermore, in CSAMT/CSRMT setups the transmitters are grounded, making their numerical treatment more difficult comparing to the conventional audio and radio magnetotellurics (AMT/RMT), respectively. For these reasons, a CSAMT/CSRMT simulation easily results in a discrete problem with millions of discrete unknowns. Beyond being very large, the system matrix is extremely ill-conditioned. Thus, both preconditioned iterative solvers and sparse direct solvers are commonly considered. Although the use of unstructured tetrahedral grids [17] or hexahedral grids with hanging nodes [18] reduces the size of the discrete problem, an efficient solver still will be required.

Direct solvers are attractive for their universality and robustness [6,19], though they require very long initialization and large memory allocations. Thus, iterative solvers are more typically used in practical modeling and inversion. The commonly applied preconditioners are incomplete factorization completed by divergence correction [20] and multigrid [21,22].

In this paper, we assess performance of the two preconditioned iterative solvers introduced by the authors in [23]. They have been successfully tested on MT and marine controlled-source electromagnetic (CSEM) [23] and land CSEM [24] setups. Here, we applied the two preconditioners to quasi-stationary simulation in the CSRMT frequency range. Both serial and parallel modeling performance were evaluated. The solvers were applied to the secondary field formulation to avoid excessive gridding near the source.

The paper is organized as follows. We review forward modelling and briefly discuss implementation in Section 2. In Section 3, we verify modeling accuracy of our code versus an integral equation code [25]. Next, we present results of numerical simulation of a study area in Moscow syneclise. We compare performance of the two preconditioners in a single-threaded mode. Then we present a complete simulation over a realistic geologic model conducted on a cluster. Final remarks are given in Section 4.

2. Three-Dimensional Simulation of CSAMT/CSAMT Data

2.1. Efficient Finite-Difference Simulation

In CSAMT/CSRMT, the five components of the electromagnetic filed are measured (E_x, E_y, H_x, H_y, H_z) 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],

$$\begin{gathered} Z_{xy}=\frac{E_{x1}{H}_{x2}-E_{x2}{H}_{x1}}{H_{x1}{H}_{y2}-H_{x2}{H}_{y1}}, \\ {Z}_{yx}=\frac{E_{y1}{H}_{y2}-E_{y2}{H}_{y1}}{H_{x1}{H}_{y2}-H_{x2}{H}_{y1}}, \end{gathered}$$

(1)

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 $\mathit{E}$ and $\mathit{H}$ at several locations and frequencies.

We consider electromagnetic modeling in 3D heterogeneous isotropic media in the frequency domain. Assuming time dependence of $e^{-i\omega t}$, the electric field, $\mathit{E}\left(x,y,z\right)$, satisfies the following system of partial differential equations,

$$curl curl \mathit{E}-i\omega {\mu }_0\sigma \mathit{E}=i\omega {\mu }_0\mathit{F},$$

(2)

where $\omega$ is the source angular frequency, ${\mu }_0$ is the magnetic permeability of the free space, $\sigma \left(x,y,z\right)$ is the conductivity model, and $\mathit{F}\left(x,y,z\right)$ is the source current density. This is a linear system of three scalar equations with respect to three scalar entries of $\mathit{E}\left(x,y,z\right)$. We solve this partial-differential equation system in a bounded domain, completed by zero Dirichlet boundary conditions. The magnetic field $\mathit{H}$ is then computed via Faraday’s law,

$$\mathit{H}=\frac{1}{i\omega {\mu }_0} curl \mathit{E}.$$

(3)

Since geological models predominately vary in the vertical direction, we can always split the conductivity model into background and anomalous parts,

$$\sigma \left(x,y,z\right)={\sigma }_b\left(z\right)+ {\sigma }_a\left(x,y,z\right).$$

(4)

We further assume that the following double inequality holds,

$$\alpha {\sigma }_b\left(z\right)\le \sigma \left(x,y,z\right)\le \beta {\sigma }_b\left(z\right),0<\alpha \le 1\le \beta ,$$

(5)

which controls the contrast of anomalous conductivity with respect to the background.

Given a non-uniform computational grid with $N_x\times N_y\times N_z$ cells, we apply the conventional edge-based second-order FD discretization. The FD method approximates the unknown electric field $\mathit{E}=\left(E_x,E_y,E_z\right)$ with a finite set of discrete values $\left\{E_{i+\frac{1}{2} j k},E_{i j+\frac{1}{2} k},E_{i j k+\frac{1}{2}}\right\}$ [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 $\mathit{e}$ of the discrete electric fields $\left\{E_{i+\frac{1}{2} j k},E_{i j+\frac{1}{2} k},E_{i j k+\frac{1}{2}}\right\}$ introduced above. Now we can write the FD discretization of (2) in a matrix form:

$$\mathit{A} \mathit{e}=i\omega {\mu }_0\mathit{f},$$

(6)

where $\mathit{A}$ 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 $n$, $n\approx 3N_x{N}_y{N}_z$.

Let ${\mathit{A}}_b$ 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):

$${\mathit{A}}_b^{-1}\mathit{A} \mathit{e}=i\omega {\mu }_0{\mathit{A}}_b^{-1}\mathit{f}.$$

(7)

We will refer ${\mathit{A}}_b^{-1}$ as the Green’s function (GF) preconditioner. The complexity of applying the GF preconditioner is $O\left(n^{4/3}\right),$ and auxiliary memory required is near $3n$ only. Applying the analysis presented in [23] the condition number of the GF preconditioned system can be estimated as follows,

$$cond\left({\mathit{A}}_b^{-1}\mathit{A}\right)\le \raisebox{1ex}{$\beta $}\!\left/ \!\raisebox{-1ex}{$\alpha $}\right..$$

(8)

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 $\mathit{\Sigma }$ and ${\mathit{\Sigma }}_b$ 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:

$${\mathit{𝓖}}_b^M=2i\omega {\mu }_0{\mathit{\Sigma }}_b^{\frac{1}{2}}{\mathit{A}}_b^{-1}{\mathit{\Sigma }}_b^{\frac{1}{2}}+\mathit{I},$$

(9)

$${\mathit{K}}_1=\frac{1}{2}\left(\mathit{\Sigma }+{\mathit{\Sigma }}_b\right){\mathit{\Sigma }}_b^{-\frac{1}{2}},\hspace{1em}{\mathit{K}}_2=\frac{1}{2}\left(\mathit{\Sigma }-{\mathit{\Sigma }}_b\right){\mathit{\Sigma }}_b^{-\frac{1}{2}}.$$

(10)

Using this operator, Equation (7) can be written in an equivalent form as follows (see Appendix A):

$$\hat{\mathit{e}}={\mathit{𝓖}}_b^M {\mathit{K}}_2{\mathit{K}}_1^{-1}\hat{\mathit{e}}+i\omega {\mu }_0{\mathit{\Sigma }}_b^{\frac{1}{2}}{\mathit{A}}_b^{-1}\mathit{f},$$

(11)

where $\hat{\mathit{e}}={\mathit{K}}_1\mathit{e}$. By introducing a new operator, $\mathit{C}={\mathit{𝓖}}_b^M {\mathit{K}}_2{\mathit{K}}_1^{-1}$, we rewrite Equation (11) as follows:

$$\left(\mathit{I}-\mathit{C}\right)\hat{\mathit{e}}=i\omega {\mu }_0{\Sigma }_b^{\frac{1}{2}}{\mathit{A}}_b^{-1}\mathit{f}.$$

(12)

We will refer this system a contraction operator (CO) preconditioned system, since it was shown in [23] that $\Vert \mathit{C}\Vert <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,

$$cond\left(\mathit{I}-\mathit{C}\right)\le \max \left\{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\alpha $}\right.,\beta \right\}.$$

(13)

The complexity of applying the CO preconditioner is $O\left(n^{4/3}\right)$ 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 $\mathit{E}\left(x,y,z\right)$ is represented as a sum ${\mathit{E}}_a\left(x,y,z\right)+{\mathit{E}}_b\left(x,y,z\right)$, where ${\mathit{E}}_b\left(x,y,z\right)$ is the response due to background conductivity model known analytically [30,31]. In this case, the secondary field ${\mathit{E}}_a\left(x,y,z\right)$ will be the unknown function and its FD discretization is performed. The actual source $\mathit{F}\left(x,y,z\right)$ in Equation (2) is substituted with the secondary source ${\sigma }_a\left(x,y,z\right){\mathit{E}}_b\left(x,y,z\right)$.

2.2. 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 $i\omega {\mu }_0\sigma$ term in Equation (1) with $i\omega {\mu }_0\sigma -{\omega }^2{\mu }_0{\epsilon }_0{\epsilon }^{rel}$, where ${\epsilon }_0$ is the vacuum dielectric permittivity and ${\epsilon }^{rel}\ge 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).

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. The background 1D resistivity model.

#	Top, m	Thickness, m	Resistivity, Ωm
1	$-\infty$	$\infty$	108
2	0	8	500
3	8	92	20
4	100	10	104
5	110	10	20
6	120	$\infty$	104

. Computations were conducted in the quasi-stationary regime, where the squared wavenumber of $i$-th layer was given by $i\omega {\mu }_0{\sigma }_i$. The air conductivity was set to ${\sigma }_0$ = ${10}^{-8}$ S/m. In the second case, the squared wavenumbers of each $i$-th layer were defined as $i\omega {\mu }_0{\sigma }_i-{\omega }^2{\mu }_0{\epsilon }_0{\epsilon }_i^{rel}$. We set ${\epsilon }_0^{rel}=1$ in the air, and ${\epsilon }_i^{rel}=4, i=1,2,\dots$ 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.

2.3. 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.

3. Numerical Experiments

3.1. 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 parallelepiped body had coordinates of the two opposite corners (110, −28, 20) m and (240, 28, 52) m; thus, that was a brick 132 × 56 × 32 m³ with its top at 20 m (Figure 3).

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³. It included both the anomalous body and receiver location (Figure 4).

The grid step size was identical in all frequencies, h = 2 m. The padding part of the grids was different for different frequencies. The overall physical size of the grids varied from 56.0 × 56.0 × 28.1 km³ at 0.1 kHz to 6.7 × 6.4 × 3.3 km³ at 10 kHz. The number of the internal edges was between 17 M at 10 kHz and 23 M at 0.1 kHz.

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.

We observed a close match between individual components. The normalized misfit between IE and FD responses, that is $\Vert u_{IE}-u_{FD}\Vert /\Vert u_{IE}\Vert$, 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.

3.2. 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.

The receiver grid consisted of 192 receiver stations and two orthogonal transmitter lines (see Figure 6 and Figure 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 top part of the column (above 93 m) is composed of mostly clayey sediments with a resistivity of approximately 20 Ωm except for the shallowest Quaternary sands. The sands have a resistivity of about 500 Ωm. There are two thin (~10 m) layers of limestones and dolomites at the depth of 93 m and 115 m. Their resistivity is estimated to 5000 Ωm. The depth of the top of the first carbonate layer was clearly detected by CSRMT across the area. It is underlain by Famennian rocks that consist mostly of dolomites and marls with the overall thickness of 140 m. This formation has a resistivity of about 10 Ωm. Since this estimate came from AMT data and the Famennian formation contains several thin dolomite layers known to be resistive within the region, this resistivity estimate can be regarded as the lateral resistivity. Below this layer, the well entered hard limestones showing its upper 30 m. We estimated the thickness of this layer as 30 m and resistivity as 1000 Ωm, based on the results of the AMT inversion and numerical simulation experiments. Properties of underlying rocks down to 1 km have been obtained from 1D inversion of a single AMT curve. The depth of the Famennian rocks was estimated as 480 m. A very conductive layer (1.5 Ωm) in the range from 480 to 730 m, as well as a resistive (630 Ωm) basement below 730 m, were introduced by the 1D inversion. The final 1D reference model is presented in

Table 2. The reference 1D resistivity model.

#	Top, m	Thickness, m	Resistivity, Ωm
1	$-\infty$	$\infty$	108
2	0	12	18
3	12	32	20
4	44	26	10
5	70	23	16
6	93	10	5000
7	103	12	11
8	115	17	5000
9	132	140	11
10	272	32	1000
11	304	176	11
12	480	250	1.5
13	730	$\infty$	670

.

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.

3.3. 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. Parameters of the finite-difference grids.

Grid Parameters	192 Hz	320 Hz	576 Hz
FD grid dimensions	114 × 110 × 151	136 × 128 × 183	166 × 158 × 138
Num of discrete unknowns	5.6 M	9.4 M	10.7 M
Grid step size in core domain, m	7.4	5.7	4.3
Core domain, m	526 × 497 × 807	534 × 488 × 800	527 × 492 × 400

.

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⁻¹⁰. Performance of the two preconditioned iterative solvers is presented in

Table 4. Comparison of Green’s function (GF) and contraction operator (CO) preconditioned iterative solvers.

Solver	Solver Step	192 Hz	320 Hz	576 Hz
CO	RHS computation, sec	2994	5762	9239
	Iteration count	346	389	480
	Iterative solver, sec	3658	8177	13,673
	Single iteration, sec	10.6	21.0	28.5
GF	RHS computation, sec	2993	5759	9258
	Iteration count	3042	3713	5000 (*)
	Iterative solver, sec	31,932	77,253	141,133 (*)
	Single iteration, sec	10.5	20.8	28.2

(*) Did not converge in 5000 iterations.

. 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 air–ground 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. Computational grids at different frequencies.

Frequency, Hz	Core Domain, m	Grid Step in Core Domain, m	Discrete Unknowns
192	526 × 497 × 800	3.00	11.7 M
320	534 × 488 × 800	3.00	16.2 M
576	527 × 492 × 397	2.99	14.0 M
960	524 × 497 × 300	3.00	10.9 M
1500	525 × 493 × 200	2.99	11.0 M
2500	530 × 489 × 150	3.00	13.5 M
3500	531 × 490 × 150	3.00	17.1 M
5500	529 × 490 × 147	2.73	25.2 M
7500	529 × 491 × 148	2.34	34.5 M
15,000	528 × 492 × 148	1.67	70.2 M
25,000	528 × 493 × 150	1.95	47.5 M
35,000	530 × 491 × 147	3.26	17.7 M
55,000	528 × 491 × 97	2.56	23.4 M
75,000	529 × 493 × 78	2.22	29.2 M
150,000	530 × 492 × 38	1.54	45.9 M
250,000	529 × 492 × 40	1.21	76.0 M
350,000	528 × 493 × 19	1.00	86.9 M
550,000	529 × 492 × 10	0.91	93.0 M

.

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. Performance of the forward simulation with OpenMP parallelization on 24 cores.

Frequency, Hz	Discrete Unknowns	RHS Computation, Sec	Iteration Count	Iterative Solver, Sec
192	11.7 M	398	459	732
320	16.2 M	572	479	1076
576	14.0 M	729	460	1043
960	10.9 M	318	502	845
1500	11.0 M	380	535	934
2500	13.5 M	561	443	1052
3500	17.1 M	765	447	1568
5500	25.2 M	1208	665	3555
7500	34.5 M	1795	553	4700
15,000	70.2 M	4489	578	12,741
25,000	47.5 M	3006	619	7264
35,000	17.7 M	768	495	1766
55,000	23.4 M	1122	436	2269
75,000	29.2 M	1362	227	1947
150,000	45.9 M	2304	64	956
250,000	76.0 M	4140	21	855
350,000	86.9 M	4120	17	1048
550,000	93.0 M	3994	14	1108

.

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.

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).

4. 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.