HPC Geophysical Electromagnetics: A Synthetic VTI Model with Complex Bathymetry

Abstract

We introduce a new synthetic marine model for 3D controlled-source electromagnetic method (CSEM) surveys. The proposed model includes relevant features for the electromagnetic geophysical community such as large conductivity contrast with vertical transverse isotropy and a complex bathymetry profile. In this paper, we present the experimental setup and several 3D CSEM simulations in the presence of a resistivity unit denoting a hydrocarbon reservoir. We employ a parallel and high-order vector finite element routine to perform the CSEM simulations. By using tailored meshes, several scenarios are simulated to assess the influence of the reservoir unit presence on the electromagnetic responses. Our numerical assessment confirms that resistivity unit strongly influences the amplitude and phase of the electromagnetic measurements. We investigate the code performance for the solution of fundamental frequencies on high-performance computing architectures. Here, excellent performance ratios are obtained. Our benchmark model and its modeling results are developed under an open-source scheme that promotes easy access to data and reproducible solutions.

1. Introduction

The Earth’s subsurface holds natural resources which are fundamental for local and regional development. Obtaining accurate images of water reservoirs, mineral, and energy sources deep below the surface is a crutial step for their management and exploitation. Geophysical imaging allows us to obtain detailed maps of the Earth’s interior. This is achieved by analyzing the deformations and electromagnetic (EM) fields measured at the surface. The EM modeling routines estimate the EM fields arising from induced electric currents in the Earth’s subsurface. The EM response to that excitation source depends on the electrical distribution of geological properties. From this dependence, it is possible to extract useful subsurface information to improve and reinforce reservoirs’ characterization and interpretation [1]. In consequence, EM approaches are now a well-established tool in geophysics, with applications in a wide range of fields including hydrocarbon and mineral exploration, reservoir monitoring, CO${}_2$ storage characterisation, geothermal reservoir imaging, water prospecting, and more.

Geophysical EM modeling is an active and dynamic research area. The most often used methods for resolving Maxwell’s equations are the Integral Equation method (IE; [2,3,4,5]), and different versions of the differential Equation (DE) method, such as Finite Difference (FD; [6,7,8,9,10,11]), Finite Elements (FE; [12,13,14,15]), and Finite Volume (FV; [16,17,18]). Furthermore, several approaches to accelerate the solution of Maxwell’s equations have been evaluated in different application contexts. Out of these approaches, the conjugate gradient method [19], vectorization techniques [20], multifrontal methods [13], and sparse matrix reordering schemes [21], stand out.

Several 3D modelers have been developed for mapping the subsurface through the study of the electrical conductivity/resistivity as a diagnostic physical property. Out of these modeling tools, custEM [22], emg3d [23], PETGEM [24], SimPEG [25], stand out. These modeling routines should be particularly sought for: (a) providing accurate solutions in a feasible runtime; (b) tackling problems efficiently; (c) bringing flexibility to cope with a variety of real-life models. The improvement of these modeling tools have a critical role in solving the next generation of geoscience challenges. These problems are complex, multidisciplinary, and require collaboration to understand and solve the physical equations, pre-process and post-process the associated data with physical experiments, and build interpretations from the analysis of the numerical results. These are the principal motivations of our work, together with the necessity for more reproducible modeling results.

This paper introduces a new model for 3D controlled-source electromagnetic method (CSEM) surveys. The proposed model includes large conductivity contrast with anisotropic (vertical transverse isotropy, VTI), and a complex bathymetry profile, making it challenging and ideal for addressing the topic of validation and for testing new algorithms and modeling routines. We have used the PETGEM code to compute the synthetic EM responses, which has proven to be an efficient large-scale modeling routine on cutting-edge high-performance computing (HPC) clusters. We acknowledge that the CSEM modeling methodology presented in this paper is well-established. However, the availability of data (resistivity model, input mesh, electromagnetic responses) for reproductively/verification purposes continues to be a limiting issue in the geophysical EM modeling community. To reverse this situation, the proposed model and numerical results are based on an open approach. Furthermore, we state that our open data and benchmark are valid and useful for different application contexts (e.g., oil and gas, geothermal exploration, CO${}_2$ sequestration, among others) and passive-source EM schemes such as the Magnetotelluric method (MT). This transversality feature favors the expertise transfer from mature energy applications (e.g., oil and gas) to emerging energy applications (e.g., geothermal exploration).

The rest of the paper is organized as follows. Section 2 provides a comprehensive description of the mathematical background of the CSEM forward modeling. In Section 3, we present details about the PETGEM code and its HPC workflow for EM modeling. In Section 4, we provide a detailed description of the proposed marine CSEM model. In addition, we perform several PETGEM simulations to investigate the EM field pattern and the code performance on HPC architectures. We conclude with a discussion in Section 5. We provide summary remarks in Section 6.

2. EM Modeling Theory

The EM forward modeling problem is mathematically described by the frequency-domain Maxwell’s equations in diffusive form, written as [26]

$$\begin{array}{cc}\hfill \nabla \times \mathbf{E}& =i\omega \mu \mathbf{H}+\mathbf{K},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \nabla \times \mathbf{H}& =\mathbf{J}+(\sigma +i\omega ϵ)\mathbf{E},\hfill \end{array}$$

(2)

where $\mathbf{E}$, $\mathbf{H}$ are the electric and magnetic fields, respectively; $\mathbf{J}$, $\mathbf{K}$ are the electric and magnetic sources, respectively; i is the imaginary unit; $\omega$ denotes the angular frequency; $\mu$ is the constant free-space magnetic permeability due Earth materials are not magnetizable; $ϵ$ denotes the constant model permittivity; and $\sigma$ is the variable electric conductivity tensor. According to the type of the excitation source, the EM forward problem may be classified as active-source (e.g., CSEM) and passive-source (e.g., MT) methods. In this paper, we focus on the CSEM method and we provide a brief outline and details essential for understanding it. For a comprehensive introduction to the CSEM and MT methods, we refer to [24,27,28,29].

In the field of CSEM modeling, it is commonly to formulate the Equations (1) and (2) in terms of total electric field. After substituting Equation (1) into Equation (2), we obtain

$$\begin{array}{c}\hfill \nabla \times \nabla \times \mathbf{E}-i\omega \mu \sigma \mathbf{E}=i\omega \mu \mathbf{J},\end{array}$$

(3)

which is known as the curl-curl formulation of the problem in terms of the total electric field. In Equation (3) we impose the usual assumption $\sigma \gg \omega ϵ$. When CSEM methods are considered, the electric source $\mathbf{J}$ corresponds to a electric dipole oscillating at a single frequency. Furthermore, no magnetic sources are generated, thus $\mathbf{K}=\mathbf{0}$.

The numerical approximation of the EM fields by the High-order Edge Finite Element Method (HEFEM) requires the variational form of Equation (3). For more details and proofs about this high-order numerical discretization, we refer to [24,30,31]. For our modeling purposes the computational domain is truncated by homogeneous Dirichlet boundary conditions.

3. HPC Workflow for EM Modeling

We have implemented the above algorithm and numerical schemes in a fully distributed fashion using Python 3 language, mpi4py, and petsc4py packages. The net result is the PETGEM code, which is an EM tool to solve both active-source and passive-source EM methods in 3D arbitrary marine/land problems under anisotropic conductivities. A triple helix model based on global polynomial variants (p-refinement), tailored tetrahedral meshes (h-refinement), and massively parallel computations, stand out as the main distinguishing factors of PETGEM with respect to other modeling routines.

The overall PETGEM software stack is shown in Figure 1. To point out the virtues of this code environment, we summarize the essential aspects its main modules:

• Provided by user: this set of modules provides functionalities to parse user parameters regarding the physics of the EM problem to be solved (e.g., frequency and conductivity model), the solver type (e.g., direct or iterative), and the mesh file.
• Specific modules: depending on the source type (active or passive), these modules provide support to define the modeling work-flow. They are in charge of data preprocessing (e.g., import mesh file and compute its associated data structures such as dof connectivity) and data postprocessing (e.g., computation of EM responses on a set of points and output of files for posterior analysis).
• Generic modules: Once the modeling’s specifics are provided, these modules perform the assembly and solution of the linear system of Equations (LSE). These tasks are entire parallel by using the mpi4py and petsc4py packages.

The software stack mentioned above is modular, simple, and flexible. Furthermore, its HPC support make it appropriate to be used in an inversion procedure scheme (e.g., iterative or stochastic). For clarity, Figure 2 depicts the block diagram overview of the PETGEM code. This diagram can be summarize as follows:

• Following the input parameters, a set of data are preprocessed (mesh, conductivity model and receivers positions).
• A problem instance is created.
• The domain decomposition is performed, and the main data structures are created.
• Parallel assembling of the LSE ($Ax=b$).
• The LSE is solved in parallel by calling a ksp PETSc object.
• Interpolation of electromagnetic responses and post-processing parallel stage.

We state that our modeling work-flow is open-source under the BSD-3 license. Please notice that the majority of the current tools for 3D geo-electromagnetic modeling lack this last feature, which makes it hard for other researchers to assess whether their underlying formulation and algorithmics suit whichever computing architecture those tools are deployed on. However, our parallel code is entirely open-source. Thus, it allows for a complete exploration of its numerical formulation and implementation. This feature eases its efficient deployment on different computational architectures. PETGEM code is hosted and versioned in a public on-line repository with unit-testing and continuous integration control. This open model engages users to submit and track issues, promoting the building of global communities around not only PETGEM project but also for a broader audience since it fosters research at the intersection of EM modeling, high-order numerical methods, and HPC.

4. Numerical Experiments

We design a new synthetic 3D CSEM model that integrates relevant features for the EM community. More concretely, the model under consideration includes large conductivity contrast with VTI, complex seabed profile, and realistic domain dimensions. There are no semi-analytical solutions for such a resistivity model to verify our synthetic EM responses, and we can only corroborate them by comparing different numerical approximations. On the off chance that distinctive discretizations and implementations of Maxwell’s equations surrender the same result, it gives certainty in their numerical precision. Then, for all experiments, we consider a normalized root-mean square difference (NRMSD) between two synthetic responses $Q_1$ and $Q_2$, given by

$$\begin{array}{c}\hfill \mathrm{NRMSD}(\%)=\frac{|Q_1-Q_2|}{\left(\right|Q_1|+|Q_2\left|\right)/2},\end{array}$$

(4)

to which we allude as basically the normalized difference through the paper.

The simulations have been performed on Marenostrum supercomputer, which is based on Intel Xeon Platinum processors from the Skylake generation. It is a Lenovo system composed of SD530 Compute Racks, an Intel Omni-Path high performance network interconnect and running SuSE Linux Enterprise Server as operating system. This general purpose supercomputer consists of 48 racks housing 3456 nodes with a grand total of 165,888 processor cores and 390 Terabytes of main memory. Compute nodes are equipped with:

• Two socket Intel Xeon Platinum 8160 CPU with 24 cores each at 2.10 GHz for a total of 48 cores per node
• L1d 32 K; L1i cache 32 K; L2 cache 1024 K; L3 cache 33,792 K
• Total of 96 GB of main memory 1.880 GB/core, 12 × 8 GB 2667 Mhz DIMM (216 nodes high memory, 10,368 cores with 7.928 GB/core)
• A 100 Gbit/s Intel Omni-Path HFI Silicon 100 Series PCI-E adapter
• A 10 Gbit Ethernet
• A 200 GB local SSD available as temporary storage during jobs

The processors support well-known vectorization instructions such as SSE, AVX up to AVX–512. We point out that in the corresponding PETGEM website, the data model (e.g., geometry, mesh file, parameters file) and modeling results can be download for comparison purposes.

4.1. 3D VTI Marine Model with Bathymetry (3D-VTI-B)

We propose a synthetic 3D VTI model that includes large conductivity contrast and challenging bathymetry profile. The model, depicted in Figure 3, is comprised of an air layer (${\rho }_{\mathrm{air}}=1^9\mathsf{\Omega }$m) with constant thickness of $2.2$km. The seawater is an isotropic medium (${\rho }_{\mathrm{water}}=0.3\mathsf{\Omega }$m) where the sea-bottom depth ranges from $289$m to $2617$m. The subsurface consists of one sedimentary formation with VTI (${\rho }_{\mathrm{seds}\left(h\right)}=0.5\mathsf{\Omega }$m and ${\rho }_{\mathrm{seds}\left(v\right)}=0.25\mathsf{\Omega }$m). Finally, a resistivity unit (${\rho }_{\mathrm{reservoir}\left(h\right)}=1000\mathsf{\Omega }$m and ${\rho }_{\mathrm{reservoir}\left(v\right)}=200\mathsf{\Omega }$m) denoting hydrocarbon is embedded at the sediments layer. The horizontal extent of the hydrocarbon reservoir is $1.2$km × $2.5$km. The model encompasses a $22.8$km × $17.4$km × $7.5$km volume with its origin at $x=718.6$km, $y=1152.6$km, $z=-2.5$km. We use an x-oriented electric dipole at $2$Hz located above the seafloor ($x=728.5$km, $y=1161$km, $z=0.975$km). The transmitter moment is $1$Am.

We based our computational model design on the skin-depth ($\delta$) principle, defined as the effective depth of penetration of EM energy in a conducting medium, where the amplitude of a plane wave in a whole space has been attenuated to $1/e$ or $37\%$ [24]. For this modeling test, the $\delta$ value is approximately $11$km in terms of the host resistivity of reservoir. By using this parameter, and following the rules proposed by [31], seventy-five measuring sites are located along the x-axis in the range $x=[725,735]$km with a receiver distance of $135$m, directly above the seafloor ($y=1152.6$ km). Given the resistivity model and reservoir dimensions, we state that this measuring sites configuration satisfies the accuracy and detection requirements of the resistive body. Furthermore, we design an $hp$-adapted mesh for basis order $p=2$. The resulting mesh consists of 1,918,108 degrees of freedom (dof).

4.2. EM Fields Analysis

To investigate the impact of the reservoir presence on the measured EM responses, we compare the numerical solutions for a homogeneous model (${\rho }_{\mathrm{seds}}={\rho }_{\mathrm{reservoir}}$) against a non-homogeneous model (${\rho }_{\mathrm{seds}}\ne {\rho }_{\mathrm{reservoir}}$, which is equivalent to the initial configuration). We compute custEM and PETGEM solutions. A multifrontal parallel solver MUMPS has been used to solve the resulting LSE for the unknown EM fields. The MUMPS solver is supported by both EM modeling routines.

For both homogeneous and non-homogeneous models, Figure 4 shows the obtained electric field amplitude $|{\mathbf{E}}_x|$ and phase ${\varphi }_x$ along the in-line receiver profile. Here, it can be seen that both numerical approximations are similar for receivers located far from the resistivity unit vicinity ($x=725$km to $x=730$km). However, for the case where the reservoir is present, the electric field amplitudes and phases are only modified close the resistivity unit, which corresponds to a local effect ($x=730.5$km to $x=735$km). More concretely, for the non-homogeneous case, the $|{\mathbf{E}}_x|$ and ${\varphi }_x$ are modified by different orders of magnitude along with the in-line profile. An excellent agreement between the custEM and PETGEM solutions can be seen. The synthetic EM responses have mostly a relative error of less than 1–2%. The cross-validation of the numerical approximations for non-homogeneous yields an excellent agreement with NRMSD within a few per cents. For clarity, Figure 5 shows the obtained NRMSD for amplitude and phase responses. This EM pattern is as expected and has been broadly studied in previous works for several and different CSEM setups [28,29]. Nonetheless, we acknowledge that the electric field pattern depends on the input model setup (e.g., source type, fundamental frequency, resistivity, depth of station where EM fields are measured). For a comprehensive analysis of the impact of these factors on EM field behaviour we refer to [32].

We compute the electric field amplitude around the reservoir location to present a full analysis of the 3D EM behavior. Figure 6 shows the EM responses for three view profiles. The amplitude distributions depicted in the first column of Figure 6 illustrate filed perturbations due to the presence of the resistivity unit. The second column of Figure 6 shows additional electric field perturbations due to the VTI. Finally, the third column of Figure 6 shows the amplitude variations caused by the reservoir presence. We point out that the significant differences in the electric field components imply that the vertical EM field is more susceptible to the perturbations of the resistivity model and is more informative than the in-line EM field for laterally heterogeneous targets. Then, this feature could be employed to get better signal-to-noise trade-offs on stations near to the vicinity of the resistivity unit [32]. Nonetheless, we acknowledge that this EM field pattern is reduced and only effective close to the reservoir location ($x=730.5$km to $x=735$km). Thus, its effect depends on where the reservoir unit is located with respect to the source. Nevertheless, the effect on EM measurements should be considered to avoid misunderstanding interpretations [32].

4.3. Performance Analysis

As a second experiment, we study the parallel performance of the PETGEM code for the solution of the proposed 3D-VTI-B model. For this test, we solve five fundamental frequencies on adapted meshes for basis order $p=2$. The obtained grid statistics are summarized in

Table 1. Statistics of adapted meshes and performance results for 3D-VTI-B model. The frequency (Hz), number of elements, number of dof, run-time (minutes), speed-up S, and parallel efficiency E (%) for five relevant simulations are given.

Frequency	Elements	Dof	Run-Time			S			E
			48	528	1008	48	528	1008	48	528	1008
10	602,593	3,678,324	88.85	9.06	4.99	-	9.8	17.79	-	89.09	84.76
5.0	615,378	3,737,031	76.85	8.35	4.54	-	9.2	16.90	-	83.63	80.47
2.5	611,482	3,714,628	68.37	7.51	4.19	-	9.1	16.29	-	82.72	77.61
1.0	613,247	3,721,364	61.45	7.02	4.12	-	8.74	14.89	-	79.54	70.95
0.5	608,577	3,864,769	49.78	6.10	3.61	-	8.15	13.77	-	74.09	65.61

. We state that the design of adapted meshes can be advantageous for the solution of the CSEM test under study. A close assessment of Table 1 shows that the number of tetrahedal elements and dof in the mesh remains constant for all frequencies. Further, the run-time to reach the solution for each fundamental frequency is almost constant. The iterative GMRES solver provided by PETSc has been used to solve the LSE.

We investigate two main performance metrics on distributed-memory clusters, namely the speed-up ratios and parallel efficiency of the code. We analyze how the run-time varies with the number of CPU N (e.g., 48, 528, 1008 CPUs) for a fixed total problem size. For more details about these performance metrics, we refer to [31]. The obtained performance ratios are shown in Table 1. Figure 7 shows the performance indicators for each target frequency. We achieved a near linear speed-up ratio for up to 528 processors. For a higher number of processors, the performance decreases mainly because the execution becomes dominated by the communication task during the solver phase. However, the scalability ratio keeps increasing up to 1008 CPUs.

Furthermore, in our numerical results, better performance metrics for high-frequency solutions can be observed. We point out that these performance ratios depend primarily on the parallel solver features. Because we employ PETSc solver implementations, no particular work was undertaken to reduce run-time, as this is an completely different goal. However, these numerical results indicate a remarkable benefit when parallel executions are employed (similar conclusions to that described in [31,32]).

5. Discussion

The development and implementation of routines for 3D geo-electromagnetic modeling has expanded within the last decade. Consequently, nowadays, there are diverse tools available to solve arbitrarily models of the 3D geo-electromagnetic modeling. However, most of the existing routines for computation of synthetic EM responses absence an open-source environment, which makes it difficult for other EM modellers or developers to study, adjust, and expand the code capabilities to their own requirements. Furthermore, most of the data-sets are proprietary or not available in open repositories for free use. Previous ideas are the principal motivations for this paper and the introduction of a 3D-VTI-B model.

The proposed model is realistic in terms of dimensions and physical parameters (e.g., resistivity model with high contrasts and large range of fundamental frequencies). We compute the solutions for the 3D-VTI-B model. Overall, we obtain an excellent match between the numerical approximations obtained with custEM and PETGEM. The main purpose of the first test is to study the impact of the reservoir presence on the measured EM responses. The cross-validation between custEM and PETGEM yields a similar EM pattern (the obtained synthetic EM responses have mostly a relative misfit of less than 1–3%). Thus, we consider these numerical results correct because comparing different modeling routines that use different numerical schemes is proper to address the topic of verification.

The 3D-VTI-B model simulations performed on tailored meshes yields a positive impact in terms of number of dofs required to satisfy a quality criterion (e.g., threshold error in synthetic EM responses). In compliance with one of the principal outcome of this paper, the grid design is a difficult and laborious task, but its expertise is critical for time to come modeling rouyine developments (e.g., implementation of 3D EM inversion algorithms). We design tailored meshes based on the skin-depth parameter as the principal quality factor to decide the characteristic grid sizes. Such griding strategy consider the behavior of EM fields computation with its diffusive pattern. The main findings of this work corroborate that tailored meshes are required to improve the algorithm’s adaptability and supply accurate EM responses in a viable computation time. The numerical results of the 3D-VTI-B model corroborate that the strict employment of tailored griding rules can be effective for 3D geo-electromagnetic modeling (e.g., the number of unknows and run-time is constant for each fundamental frequency, see Table 1). On top of that, we encourage the design and use of tailored meshes by the EM modellers.

We state that although the model we present is synthetic, it is relevant for the geo-electromagnetic community and can be used as benchmark for new modeling routines. We also point out that, given the results of the numerical tests, there is no reason to think that the proposed algorithm does not work well with real data beyond the model limitations.

Regarding the computational performance, the chosen solver is the principal discriminator. We employ the iterative GMRES solver to compute the numerical solutions for the introduced resistivity model. In general, a considerable run-time improvement can be seen by performing a strong scaling test. We state that the improvement of reported run-time has been achieved through HPC, a differentiating PETGEM feature concerning the rest of EM modeling tools. Our performance study shows that high-frequencies achieve better performance ratios than low-frequencies.

6. Conclusions

We have presented a new 3D VTI model (3D-VTI-B) that includes significant conductivity contrast and a challenging bathymetry profile. We, carried out several simulations to analyze the EM field behavior and the code performance. The relevance of these studies is based on two main arguments. First, the solution to the next generation of geoscience problems requires accurate, flexible, and efficient modeling routines. Thus, to address the validation of the computational tools, models that exhibit realistic setups are needed. Second, it is well known the scarcity of more open and reproducible numerical approximations in the realm of geophysical EM. The challenges mentioned above are the principal motivations of this work.

Because the proposed model has no semi-analytical solution, we validate our numerical approximations by comparing different synthetic EM responses. Then, we compare the computed fields of two different open-source 3D EM modeling routines. We study the resistivity unit effect by considering homogeneous and non-homogeneous configurations. The cross-validation of both synthetic EM responses shows a NRMSD of a few per cent (1–2%). For an in-depth analysis, we have generated 3D views of the EM field pattern. Numerical approximations corroborate that reservoir presence strongly alters EM fields, but its impact is restricted to the near region of the resistivity unit. Furthermore, in our experiments, the vertical EM field is more sensitive to the variations of the resistivity model and, in consequence, more informative than the in-line EM responses. Still, the pattern of EM fields should be investigated and analyzed to avoid misunderstanding conclusions.

PETGEM code offers excellent computational performance ratios. The obtained performance metrics are consistent with previously published results. Furthermore, our previously published tailored gridding strategy has been proved in the presence of complex bathymetry and for modeling different fundamental frequencies. Our mesh design approach demonstrated that can deal with challenging bathymetry and realistic physical parameters.

The whole set of these numerical results presented in this work show that PETGEM code features fulfill modeling requisites of practical setups in the context of HPC geophysical electromagnetics. We acknowledge that validating and testing 3D algorithms and codes are difficult tasks. Then, it is fundamental to have benchmark models that are developed under an open-source scheme that promotes easy access to data and reproducible solutions. Therefore, in the corresponding PETGEM website and Zenodo, the dataset and all results can be download for comparison purposes. We trust that these numerical dataset may be helpful for the geophysical community interested in HPC geoelectromagnetics. We stir the community to design more models, modeling tools, and numerical experiments under an open-source approach.