Matías W Elías
email: melias@fcaglp.unlp.edu.ar

Fabio I Zyserman

Marina Rosas-Carbajal

María Constanza

María Constanza Manassero

Three-dimensional modelling of controlled source electro-magnetic surveys using non-conforming finite element methods

Keywords: Controlled source electromagnetics (CSEM), Electromagnetic theory, Numerical modelling, Numerical solutions

INTRODUCTION

The Controlled Source Electro-Magnetic method (CSEM) is widely used in geophysical exploration. This method is based on the use of artificial Electro-Magnetic (EM) sources and the principle of EM induction: EM fields are induced within the subsurface as a response to the emitted source fields. This response is measured by a set of receivers on the surface. The method is mainly sensitive to variations of the electrical conductivity in the subsurface, in particular, bodies more conductive than their surroundings in continental regions and resistive bodies in marine regions (Liu et al., 2013). The conductivity of rocks varies considerably depending on composition, temperature, the degree of fluid saturation of their pores, and the conductivity of the saturating fluid [START_REF] Archie | The Electrical Resistivity Log as an Aid in Determining Some Reservoir Characteristics[END_REF]. The differential conductivity may also be related with mineralization events or the presence of partially molten rock in volcanic systems [START_REF] Peacock | Three-dimensional electrical resistivity model of the hydrothermal system in Long Valley Caldera, California, from magnetotellurics: MAGNETOTELLURICS OF LONG VALLEY CALDERA[END_REF].

The performance and advantages of the CSEM method in hydrocarbon exploration in marine environments have been extensively demonstrated in the literature (e.g. [START_REF] Constable | An introduction to marine controlled-source electromagnetic methods for hydrocarbon exploration[END_REF]; Constable (2010); [START_REF] Fanavoll | CSEM as a tool for better exploration decisions: Case studies from the Barents Sea, Norwegian Continental shelf, Interpretation[END_REF]; [START_REF] Attias | Controlled-source electromagnetic and seismic delineation of sub-seafloor fluid flow structures in a gas hydrate province, offshore Norway[END_REF]; [START_REF] Dunham | Application of 3D marine controlledsource electromagnetic finite-element forward modeling to hydrocarbon exploration in the Flemish Pass Basin offshore Newfoundland, Canada[END_REF]). On land, CSEM has been successfully used in the exploration of mineral resources (Sheard et al., 2005;[START_REF] An | Application of the CSAMT Method for Exploring Deep Coal Mines in Fujian Province, Southeastern China[END_REF][START_REF] Yang | Three-dimentional inversion of airborne time-domain electromagnetic data with applications to a porphyry deposit[END_REF], CO 2 storage characterization [START_REF] Girard | Time-lapse CSEM monitoring of the Ketzin (Germany) CO 2 injection using 2 x MAM configuration[END_REF][START_REF] Grayver | 3-D inversion and resolution analysis of land-based CSEM data from the Ketzin CO 2 storage formation[END_REF] and monitoring of geothermal reservoirs [START_REF] Coppo | Characterization of Deep Geothermal Energy Resources in 3-D CSEM modeling using NCFE 45 Low enthalpy sedimentary basins in Belgium using Electro-Magnetic Methods -CSEM and MT results[END_REF][START_REF] Darnet | Monitoring geothermal reservoir developments with the Controlled-Source Electro-Magnetic method -A calibration study on the Reykjanes geothermal field[END_REF].

In marine CSEM surveys, relatively high frequency energy (0.1 to 10 Hz) is emitted by a deep-towed transmitter close to the seafloor, this makes the CSEM especially sensitive to thin resistive layers embedded in a more conductive medium [START_REF] Constable | Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling[END_REF]. Since, the hydrocarbon content in the rock pores produces a strong resistive anomaly, CSEM has an advantage over the seismic method in cases where the seismic response associated with changes in physical properties is relatively small [START_REF] Da Silva | A finite element multifrontal method for 3D CSEM modeling in the frequency domain[END_REF]. In addition, seismic reflection techniques are limited in distinguishing between oil and gas-charged brine reservoirs [START_REF] Dunham | Application of 3D marine controlledsource electromagnetic finite-element forward modeling to hydrocarbon exploration in the Flemish Pass Basin offshore Newfoundland, Canada[END_REF], since small amounts of gas can result in a false seismic hydrocarbon indicator. In CSEM a reservoir with gas will be electrically conductive compared to one completely saturated with oil, thus there is a clear distinction between the gas and oil scenarios.

The CSEM survey design, and the correct interpretation of the EM data in complex geological environments, requires the use of numerical forward modeling techniques. This involves solving Maxwell's equations in order to predict the EM response by establishing the problem parameters a priori, particularly the type of excitation source and the conductivity distribution of the subsurface. Forward 3-D CSEM modeling is useful mainly when lateral variations of conductivity in the subsurface are strong [START_REF] Da Silva | A finite element multifrontal method for 3D CSEM modeling in the frequency domain[END_REF]. Regarding the transmitters used, these can be treated as electric or magnetic dipole sources [START_REF] Løseth | Electromagnetic fields in planarly layered anisotropic media[END_REF], or as a combination of them. Due to the potential of CSEM, several algorithms have been developed in the past two decades. Some popular algorithms includes those of [START_REF] Zhdanov | Integral equation method for 3D modeling of electromagnetic fields in complex structures with inhomogeneous background conductivity[END_REF]; [START_REF] Bakr | Feasibility of simplified integral equation modeling of lowfrequency marine CSEM with a resistive target[END_REF]; [START_REF] Liu | A hybrid solver based on IEM and vector FEM for 3D CSEM modeling[END_REF] based on integral equations (IE); the algorithms of [START_REF] Streich | 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy[END_REF]; [START_REF] Sasaki | Useful characteristics of shallow and deep marine CSEM responses inferred from 3D finite-difference modeling[END_REF]; [START_REF] Sommer | GPU parallelization of a three dimensional Marine CSEM code[END_REF] based on finite differences (FD), the algorithms of [START_REF] Weiss | Project APhiD: A Lorenz-Gauged A-Φ Decomposition for Parallelized Computation of Ultra-Broadband Electromagnetic Induction in a Fully Heterogeneous Earth[END_REF]; [START_REF] Jahandari | A finite-volume solution to the geophysical electromagnetic forward problem using unstructured grids[END_REF] based on finite volumes (FV); and the algorithms developed by [START_REF] Schwarzbach | Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics-a marine CSEM example[END_REF][START_REF] Da Silva | A finite element multifrontal method for 3D CSEM modeling in the frequency domain[END_REF]; [START_REF] Puzyrev | A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling[END_REF]; [START_REF] Castillo-Reyes | Parallel 3-D marine controlled-source electromagnetic modelling using high-order tetrahedral Nédélec elements[END_REF] based on finite elements (FE). These algorithms implement either the electric and/or magnetic field formulation (e.g. [START_REF] Streich | 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy[END_REF][START_REF] Da Silva | A finite element multifrontal method for 3D CSEM modeling in the frequency domain[END_REF] or the coupled-potential formulation (e.g. [START_REF] Weiss | Project APhiD: A Lorenz-Gauged A-Φ Decomposition for Parallelized Computation of Ultra-Broadband Electromagnetic Induction in a Fully Heterogeneous Earth[END_REF][START_REF] Puzyrev | A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling[END_REF] to obtain the EM response.

The most common numerical methods used to solve the CSEM forward problem are those based on differential equations (i.e. FD, FV and FE). In particular, the FD method is very attractive due to its simplicity and ease of implementation. However, it is limited by the use of regular grids, whose size largely increases with complex geological structures. Contrary to FD approaches, FV and FE methods can be applied to non-uniformly structured meshes or unstructured meshes with notable advantages [START_REF] Castillo-Reyes | Parallel 3-D marine controlled-source electromagnetic modelling using high-order tetrahedral Nédélec elements[END_REF]. IE methods have been shown to be highly efficient for low complexity models, but their computational cost increases with the complexity of the model [START_REF] Streich | 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy[END_REF]. In order to overcome this computational cost, [START_REF] Liu | A hybrid solver based on IEM and vector FEM for 3D CSEM modeling[END_REF] implemented a hybrid IE-FE method, which takes advantage of both, bounded computational domain and efficiency of the IE method, and the suitability for complex geological models of the FE method.

The FE method has become increasingly popular in EM modeling due to its accuracy and versatility in representing complex geological structures [START_REF] Um | Efficient pre-conditioned iterative solution strategies for the electromagnetic diffusion in the Earth: finite-element frequency-domain approach[END_REF][START_REF] Liu | A hybrid solver based on IEM and vector FEM for 3D CSEM modeling[END_REF]. In some of its more popular versions (which employ Nedelec's elements or edge elements), the method is advantageous because it locally enforces the divergence-free condition, which is consistent with the discontinuity of the normal electric field at interfaces across different materials [START_REF] Um | An iterative finite element time-domain method for simulating three-dimensional electromagnetic diffusion in earth[END_REF][START_REF] Da Silva | A finite element multifrontal method for 3D CSEM modeling in the frequency domain[END_REF][START_REF] Liu | A hybrid solver based on IEM and vector FEM for 3D CSEM modeling[END_REF].

Within the most widely used FE formulations, the discretization of the differential equations results in a sparse system of linear equations which is commonly solved via: 1) direct or 2) iterative solvers.

The former solver involves the explicit factorization of the sparse matrix and forward and backward elimination steps to obtain the solution. It is usually chosen for its robustness to ill-conditioning [START_REF] Um | An iterative finite element time-domain method for simulating three-dimensional electromagnetic diffusion in earth[END_REF] and its ability to obtain efficient solutions to multi-source problems [START_REF] Chung | Three-dimensional modelling of controlled-source electromagnetic surveys using an edge finite-element method with a direct solver[END_REF].

It has been successfully applied in 3-D CSEM modeling (e.g. [START_REF] Streich | 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy[END_REF][START_REF] Schwarzbach | Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics-a marine CSEM example[END_REF][START_REF] Da Silva | A finite element multifrontal method for 3D CSEM modeling in the frequency domain[END_REF][START_REF] Chung | Three-dimensional modelling of controlled-source electromagnetic surveys using an edge finite-element method with a direct solver[END_REF] overcoming the old limitations regarding high memory requirements and computational cost. Recently, da Piedade et al. ( 2021) have compared the computational cost between nodal and vector FE using a direct solver, demonstrating that the vector FE approach is the most computationally efficient. Other advances, such as the one shown in [START_REF] Zhu | An efficient parallel algorithm for 3D magnetotelluric modeling with edge-based finite element[END_REF], are focused on speeding up computation and improving the efficiency of the direct solver with an algorithm based on distributed matrix storage and three levels of parallelism.

Iterative solvers have also been used in EM modeling, particularly in 3-D CSEM (e.g. [START_REF] Um | An iterative finite element time-domain method for simulating three-dimensional electromagnetic diffusion in earth[END_REF][START_REF] Puzyrev | A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling[END_REF][START_REF] Sommer | GPU parallelization of a three dimensional Marine CSEM code[END_REF][START_REF] Dunham | Application of 3D marine controlledsource electromagnetic finite-element forward modeling to hydrocarbon exploration in the Flemish Pass Basin offshore Newfoundland, Canada[END_REF][START_REF] Castillo-Reyes | Parallel 3-D marine controlled-source electromagnetic modelling using high-order tetrahedral Nédélec elements[END_REF][START_REF] Qiu | Solution of large-scale 3D controlledsource electromagnetic modeling problem using efficient iterative solvers[END_REF], due to their relatively low storage requirements. Furthermore, iterative solvers are much easier to implement efficiently on computers with parallel-distributed memory [START_REF] Puzyrev | A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling[END_REF]. However, poor convergence behavior and spurious solutions can arise when solving the iterative method as the frequencies approach the static limit [START_REF] Peng | An accelerated solution strategy for 3D multi-frequency controlled source electromagnetic modeling using a hybrid direct-iterative solver[END_REF]. Efforts to overcome these difficulties have led to numerous strategies such as the formulation and application of novel preconditioning techniques (e.g. [START_REF] Li | An Efficient Preconditioner for 3-D Finite Difference Modeling of the Electromagnetic Diffusion Process in the Frequency Domain[END_REF][START_REF] Qiu | Solution of large-scale 3D controlledsource electromagnetic modeling problem using efficient iterative solvers[END_REF][START_REF] Peng | An accelerated solution strategy for 3D multi-frequency controlled source electromagnetic modeling using a hybrid direct-iterative solver[END_REF], the addition of a static divergence correction to the iterative process (e.g. [START_REF] Liu | Electromagnetic divergence correction for 3d anisotropic em modeling[END_REF][START_REF] Tang | Three dimensional controlled-source electromagnetic forward modeling by edge-based finite element with a divergence correction[END_REF], the modification of the original curl-curl equations with the divergence correction equation (e.g. [START_REF] Dong | Divergence-free solutions to electromagnetic forward and adjoint problems: a regularization approach[END_REF][START_REF] Li | 3D finite difference modeling of controlledsource electromagnetic response in frequency domain based on a modified curl-curl equation[END_REF], or the formulation of the governing equations in terms of potentials (e.g. [START_REF] Weiss | Project APhiD: A Lorenz-Gauged A-Φ Decomposition for Parallelized Computation of Ultra-Broadband Electromagnetic Induction in a Fully Heterogeneous Earth[END_REF][START_REF] Ye | 3-d adaptive finite-element modeling of marine controlled-source electromagnetics with seafloor topography based on secondary potentials[END_REF][START_REF] Peng | 3-D Marine Controlled-Source Electromagnetic Modeling in Electrically Anisotropic Formations Using Scattered Scalar-Vector Potentials[END_REF].

Based on a collection of nonconforming Galerking finite element methods (NCFE) for Maxwell's equations developed by [START_REF] Douglas | A Nonconforming Mixed Finite Element Method For Maxwell's Equations[END_REF], Zyserman & Santos (2000) presented a parallel iterative domain decomposed -at the differential level-algorithm to solve the 3-D forward problem in magnetotellurics (MT). Later, [START_REF] Santos | Finite element approximation of coupled seismic and electromagnetic waves in fluid-saturated poroviscoelastic media[END_REF] employed these FE developing algorithms for Pride's equations of electroseismics [START_REF] Pride | Governing equations for the coupled electromagnetics and acoustics of porous media[END_REF]. [START_REF] Zyserman | Finite element modeling of SHTE and PSVTM electroseismics[END_REF][START_REF] Zyserman | Numerical evidence of gas hydrate detection by means of electroseismics[END_REF] implemented versions of these FE to perform numerical electroseismic studies in 2D geometries. Quite recently, Manassero (2019); [START_REF] Manassero | A Reduced Order Approach for Probabilistic Inversions of 3D magnetotelluric data I: General Formulation[END_REF][START_REF] Manassero | A Reduced Order Approach for Probabilistic Inversions of 3D Magnetotelluric Data II: Joint inversion of MT and Surface-Wave Data[END_REF] used the global version of the NCFE in the development of a novel high-dimensional 3-D MT probabilistic inversion algorithm. The capabilities of the forward solver were thoroughly tested solving standard benchmark models. The global solver was proved to be more computationally efficient than other commonly used solvers. These precedents motivate our interest in the implementation of NCFE forward solvers for 3-D CSEM.

In this work, we present two different algorithms: 1) a domain decomposed iterative solver (DDFE) and 2) a global solver (GFE) to compute the direct solution. We use a secondary field formulation for both implementations. This enables an accurate representation of the different sources (electric or magnetic) in a horizontal layer model for the primary field, which is calculated semi-analytically [START_REF] Løseth | Electromagnetic fields in planarly layered anisotropic media[END_REF]. A significant reduction of the computational domain is achieved by using the first-order absorbing boundary conditions [START_REF] Sheen | Approximation of electromagnetic fields: Part I. continuous problems[END_REF][START_REF] Manassero | A Reduced Order Approach for Probabilistic Inversions of 3D magnetotelluric data I: General Formulation[END_REF]. Both implementations, DDFE and GFE, lead to codes designed to efficiently work on parallel computers, a necessary requirement in 3-D CSEM modeling to obtain solutions in practical times. This paper is organized as follows. Section 2 introduce the differential problem in terms of the primary and secondary fields, the associated variational formulation, and numerical methods of resolution. In Section 3 we describe our implementation, which includes: i) generating the primary fields for different sources in a layered media, ii) solving the global problem using the package MUlti-frontal Massively Parallel Sparse direct Solver (MUMPS), version 5.1.2 [START_REF] Amestoy | MUMPS: A General Purpose Distributed Memory Sparse Solver[END_REF] and iii) the parallel iterative procedure to solve the problem with domain decomposition (DD). Finally, in Section 4, we test our algorithms using two different resistive disk models proposed by [START_REF] Constable | Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling[END_REF] and Streich & Becken (2011a), we carry out a scalability test, and we evaluate the sensitivity of the method using a synthetic 3-D complex conductivity structure model.

3-D CSEM modeling using NCFE 7

FORWARD FORMULATION AND APPROACH

Differential model

Maxwell's equations describe the behavior of EM fields. Considering a harmonic time dependence e iωt , at a given angular frequency ω, the diffusive Maxwell's equations are

∇ × E = -iωµ 0 H, (1a)
∇ × H -σE = J 0 , (1b)
where E and H denote the electric and magnetic fields, respectively, σ the electric conductivity, µ 0 the magnetic permeability of the vacuum (we assume it constant in all space), J 0 the current distribution of the source. As usual in CSEM, displacement currents have been neglected. Conduction currents are dominant for frequencies between 10 mHz to 20 GHz in a medium where the conductivity exceeds 0.01 S/m. They are also dominant in a medium where conductivities are between 10 -6 to 10 6 S/m for frequencies less than 1 MHz [START_REF] Everett | On the physical principles underlying electromagnetic induction[END_REF]. This dominant conduction current is satisfied for the relatively low frequency range used in CSEM and the conductivity values of the subsurface typically expected in geophysical studies.

Given that the distortion of the EM field in the vicinity of the source can lead to numerical instabilities and/or to a poor representations of it, the differential model is formulated in terms of the scattered fields [START_REF] Newman | Frequency domain modeling of airborne electromagnetic responses using staggered finite difference[END_REF][START_REF] Puzyrev | Deep learning electromagnetic inversion with convolutional neural networks[END_REF]

∇ × E s + iωµ 0 H s = 0 in Ω, (2a)
σE s -∇ × H s = -(σ -σ p )E p = -F in Ω, (2b)
where Ω ⊂ IR 3 is the 3-D domain of the problem. With the aim of minimizing the effect of the artificial borders of the domain Ω, the first-order absorbing boundary conditions are used [START_REF] Sheen | Approximation of electromagnetic fields: Part I. continuous problems[END_REF][START_REF] Manassero | A Reduced Order Approach for Probabilistic Inversions of 3D magnetotelluric data I: General Formulation[END_REF])

(1 -i)P τ aE s + ν × H s = 0 on ∂Ω ≡ Γ, (3)
where a = (σ/2ωµ 0 ) 1 2 , ν is a normal unit vector to the computational border Γ and

P τ φ = φ -ν(ν • φ) = -ν × (ν × φ).
In this way, the normal field to the border Γ is absorbed, avoiding undesirable reflections. This boundary conditions also allow to reduce the size of the computational domain, since extending the boundaries far from the source to make the fields negligible is no longer necessary.

Global finite element solution

The variational formulation of Eqs.

(2) and ( 3) is obtained following [START_REF] Douglas | A Nonconforming Mixed Finite Element Method For Maxwell's Equations[END_REF]; Zyserman & Santos (2000). We test Eq. (2b) with real vector functions ϱ(x, y, z) such that ∇ × ϱ(x, y, z) is square integrable and Eq. (2a) with square integrable vector functions ζ(x, y, z) [START_REF] Santos | Global and parallelizable domain decomposed mixed FEM for 3D electromagnetic modelling[END_REF]. Eq. ( 3) is incorporated through the use of the generalized Green's formula [START_REF] Douglas | A Nonconforming Mixed Finite Element Method For Maxwell's Equations[END_REF]. Omitting the subscripts of the secondary fields to simplify the notation, the variational formulation problem reads as follows: Find (E, H) such that

(σE, ϱ) -(H, ∇ × ϱ) + (1 -i)⟨P τ aE, P τ ϱ⟩ Γ = (F, ϱ), ∀ϱ, (4a)
iωµ 0 (H, ζ) + (∇ × E, ζ) = 0, ∀ζ. (4b)
(•, •) and ⟨•, •⟩ denote inner product in Ω and on Γ ≡ ∂Ω, respectively.

Numerical solutions of Eqs. (4a)-(4b) are obtained via the nonconforming Galerkin FE method (NCFE)

of [START_REF] Douglas | A Nonconforming Mixed Finite Element Method For Maxwell's Equations[END_REF]. In the following, we comment on general aspects of the method necessary for the understanding of the formulation in this work. We refer the reader to [START_REF] Douglas | A Nonconforming Mixed Finite Element Method For Maxwell's Equations[END_REF] and references therein for an in-depth description of the method. Within nonconforming Galerkin methods, the continuity of the tangential components of the electric field is imposed only on the center of the element's faces rather than on the complete face (as imposed in traditional conforming Galerkin methods). This simplifies the calculations compared to conforming elements. The NCFE, as well as the "edge elements", also called "Whitney elements" [START_REF] Jin | The finite element method in electromagnetics[END_REF], minimize the possibility of spurious modes [START_REF] Douglas | A Nonconforming Mixed Finite Element Method For Maxwell's Equations[END_REF]. In the case of the edge elements, it was shown that their use is not enough to eliminate them completely, but this also depends on an adequate variational formulation [START_REF] Mur | The fallacy of edge elements[END_REF].

The method requires the partition of the domain Ω in N e three-dimensional (not overlapping and not necessarily homogeneous) rectangular prisms, Ω j , j = 1, ..., N e . The infinite dimensional spaces to which the solution (E, H) of Eq. ( 4a)-(4b) belong are replaced by approximating spaces of finite elements U h and V h , respectively. Table 1 shows the shape function expanding these spaces in a reference cube K = [-1, 1] 3 . It can be seen that the space expanded by the functions ζ κ is the curl of the space expanded by the functions ϱ ι . Thus, U h is expanded by the basis function set ϱ ι , ι = 1, ..., N E dof , and V h is expanded by the set of functions ζ κ , κ = 1, .., N H dof , where N E dof and N H dof are the number of degrees of freedom (dof) for the electric and magnetic field, respectively.

In this manner, the variational problem (4) is approximated in terms of the FE spaces as:

Find (E h , H h ) ∈ U h × V h such that (σE h , ϱ) - Ne j (H h , ∇ × ϱ) Ω j + (1 -i)⟨⟨P τ aE h , P τ ϱ⟩⟩ Γ = (F, ϱ), ∀ϱ ∈ U h , (5a) iωµ 0 (H h , ζ) + Ne j (∇ × E h , ζ) Ω j = 0, ∀ζ ∈ V h , (5b)
where E h ∈ U h and H h ∈ V h are the approximate electric and magnetic field in Ω. These can be expressed in terms of the basis functions as

E h ≃ N E dof ι=1 e ι ϱ ι , (6a)
H h ≃ N H dof κ=1 m κ ζ κ , (6b)
where the superscripts ι and κ are associated with the corresponding basis functions and the coefficients e ι and m κ are the unknowns to be determined. The choice of the basis functions for the magnetic field allows the coefficients m κ to be expressed in terms of the coefficients e ι (Zyserman & Santos, 2000). This yields a linear system of equations where the unknown is a vector u composed of the coefficients e ι to approximate the electric field in the entire domain. We can express it in the following matrix form

Ku = b, (7)
where K is a large, sparse, symmetric, complex valued, and ill-conditioned matrix resulting from FE discretization, also known as stiffness matrix, and b is the force vector of the secondary fields.

Solving the global linear system (7) can be computationally demanding, because of the need of large memory resources and long computing times. However, in recent years, substantial advances have been made in the development of efficient and scalable sparse matrix factorization algorithms that allow dealing with these systems (e.g da Silva et al., 2012).

Domain decomposition solution

In this section, we present the second FE algorithm to iteratively solve Eqs. (2) using domain decomposition (DD) at the differential level.

The main benefit of DD procedures is to decompose a large-scale problem (defined in a single domain)

into several small-scale problems (defined in subdomains) which can be solved in parallel [START_REF] Xue | Finite-Element Domain Decomposition Methods for Analysis of Large-Scale Electromagnetic Problems[END_REF].

Considering that the FE partition of the domain Ω matches the domain decomposition partition (non-overlapping subdomains), we can decompose the global differential problem into smaller problems j = 1, ..., N e with (E j , H j ) such that

σE j -∇ × H j = -F j in Ω j , (8a)
iωµ 0 H j + ∇ × E j = 0 in Ω j , (8b)
(1 -i)P τ aE j + ν j × H j = 0 on ∂Ω j ≡ Γ j , (8c)
where F j = F| Ω j is the restriction of F to the domain Ω j . Note that solutions to Eqs (8) are sought in each domain Ω j , = 1, ..., N e separately. In order to obtain the solution in the entire domain Ω, consistency conditions need to be imposed on all internal boundaries between subdomains. Although natural conditions are given by the continuity of the tangential components of the electric and magnetic field at the boundaries, these are replaced by the equivalent Robin-type transmission conditions [START_REF] Douglas | A Nonconforming Mixed Finite Element Method For Maxwell's Equations[END_REF], which improve efficiency of the iterative algorithm:

(ν j × H j + β jk P τ E j ) = -(ν k × H k + β jk P τ E k ) on Γ jk ⊂ ∂Ω j , (9a)
(ν k × H k + β jk P τ E k ) = -(ν j × H j + β jk P τ E j ) on Γ kj ⊂ ∂Ω k , (9b)
where Γ jk is the shared face between adjacent domains Ω j and Ω k , and β jk is a complex parameter defined on the interfaces Γ jk with positive real part and negative imaginary part.

The algorithm then separately solves the problems in each subdomain, considering common borders between adjacent subdomains. This results in an iterative procedure for calculating the solution across the entire domain. For each subdomain, problem (8) can be posed in its variational formulation using the same shape functions as in the global case, i.e., those given in Table 1.

In order to simplify the algebraic problem, we hybridize the algorithm [START_REF] Arnold | Mixed and nonconforming finite element methods: implementation, postprocessing and error estimates[END_REF]. A hybrid approach implies relaxing the continuity constraints of the tangential component of the electric field on the faces Γ jk by dissociating the values of the magnetic field in the interior of the element from the one on its faces, and incorporating a set of Lagrange (λ) multipliers associated to the tangential component of the magnetic field evaluated at the centroid of the inter-element faces. Therefore λ ≈ (ν × H j )(c jk ) on Γ jk ; where c jk is the centroid of the interface Γ jk . Notice that these Lagrange multipliers are two dimensional vectors.

Assuming an initial condition (E h,0 j , H h,0 j , λ h,0 jk , λ h,0 kj ) the DD problem reads: Find (E h,n j , H h,n j , λ h,n jk )

3-D CSEM modeling using NCFE 11 solution of the equations

(σE h,n j , ϱ) j -(H h,n j , ∇ × ϱ) j + k ⟨⟨β jk P τ E h,n j , P τ ϱ⟩⟩ Γ jk + (1 -i)⟨⟨P τ aE h,n j , P τ ϱ⟩⟩ ∂Ω j = (F j , ϱ) j + k ⟨⟨β jk P τ E h,n-1 j -λ h,n kj , P τ ϱ⟩⟩ Γ jk , ϱ ∈ U h jk , ( 10a
)
iωµ 0 (H h,n j , ζ) j + (∇ × E h,n j , ζ) j = 0, ζ ∈ V h j , (10b)
λ h,n jk = -λ h,n-1 kj + β jk (P τ E h,n-1 k -P τ E h,n j )(ξ jk ), ξ jk ∈ Γ jk , (10c)
where U h jk and V h jk denote the restrictions of the global FE spaces (U h and V h , respectively) to the domain Ω jk . The convergence of the iterative procedure to the solution (E h , H h ) of problem ( 5) is demonstrated in [START_REF] Douglas | A Nonconforming Mixed Finite Element Method For Maxwell's Equations[END_REF].

In this case, the coefficients e ι and m κ of Eqs. ( 6) are equivalent to the coefficients e ι,n+1 j and m κ,n+1 j which are determined in the iteration n + 1. In the same way as in the global case, the coefficients m κ,n+1 j can be expressed in terms of the coefficients e ι,n+1 j . This leads to a 12 × 12 linear system of the form D j e ι,n+1 j = b n j , where e ι,n+1 j are the unknowns. The vectors b n j are recalculated at each iteration, while the coefficient matrices D j remain unchanged along the iterative process.

Thus, the iterative procedure consists of:

1 Choose initial values (e ι,0 j , λ 0 jk ) for the unknowns in all cells Ω j . 2 For all domains Ω j Solve the 12 × 12 linear system for the unknowns e ι,n+1 j Compute λ n+1 jk .

3 Evaluate convergence. If it has not been achieved, go to step 2.

Convergence is achieved when the relative error of the calculated coefficients is smaller than a certain tolerance.

Sources and primary fields

The separation into primary and secondary fields allows a better representation of the field near the source, and also the representation of 1D conductivity models allows to obtain an accurate solution without large computational cost. These models can be used for sensitivity analysis before considering a three-dimensional geometry (Streich & Becken, 2011b).

The primary EM field is calculated by solving

∇ × E p + iωµ 0 H p = 0, (11a)
∇ × H p -σ p E p = J 0 , (11b)
where the conductivity σ p = σ(z) represents a horizontal layered model and the current density J 0 determines the type of source used. The sources considered in this work are electrical and magnetic point dipoles. When considering magnetic sources it is advantageous to introduce a term J 0M in Faraday's equation (Eq. 11a), this is ∇ × E p = -iωµ 0 H p -J 0M . Since magnetic sources are introduced, the source term J 0 is referred as an electric source. A point dipole electrical source can be characterized by a current line of length l = l x x + l y ŷ + l z ẑ and current amplitude I(ω) at frequency ω [START_REF] Løseth | Electromagnetic fields in planarly layered anisotropic media[END_REF]. We represent it as

J 0 = I(ω)[l x x + l y ŷ + l z ẑ]δ(x -x s ), ( 12
)
where x is an arbitrary position in Cartesian coordinates and x s = [x s , y s , z s ] is the source position. A point magnetic dipole source, characterized by its magnetic moment m = µ 0 I(ω)a, being a a vector with normal direction to the current loop and magnitude equal to its area, by analogy to the current density described in Eq. ( 12), is represented by

J 0M = -iωµ 0 I(ω)aδ(x -x s ). (13)
The above expressions represent arbitrarily oriented sources. The three sources used here: i) horizontal electric dipole (HED), ii) horizontal magnetic dipole (HMD) and iii) vertical magnetic dipole (VMD).

In a horizontally layered medium, a transformation between the space-time and frequency-wavenumber domains is convenient to solve Maxwell's equations and obtain analytical expressions of the EM fields [START_REF] Løseth | Electromagnetic fields in planarly layered anisotropic media[END_REF]. In the transmission-reflection propagation strategy proposed by the authors, the EM field is decomposed into transversal electric (TE) and transversal magnetic (TM) modes [START_REF] Ward | Electromagnetic theory for geophysical applications[END_REF]. In both TE and TM modes, the field propagates through the interfaces upwards and downwards, thus transmission and reflection coefficients between layers are calculated. The successive reflections and transmissions are interpreted as the response of the medium to the propagation of the field, which depends on the number of layers, their thickness and conductivity, and the source and receiver's vertical position. Taking into account each type of source and the EM response of the medium, explicit expressions can be obtained for the calculation of the EM field in each case. As the expressions for the three sources are given in [START_REF] Løseth | Electromagnetic fields in planarly layered anisotropic media[END_REF] in the frequency-wavenumber domain, we apply a 2D Fourier transform on the wavenumbers to obtain the fields in the space-frequency domain. In the particular case of considering isotropic horizontal layers, this leads to simple integrals in terms of Bessel functions and the polar horizontal wavenumber in cylindrical coordinates. We compute the numerical solution of these integrals using the linear digital filter technique for the Hankel transform based on an orthonormal exponential approximation of the kernel function developed by [START_REF] Guptasarma | New digital linear filters for Hankel J 0 and J 1 transforms[END_REF].

3-D CSEM modeling using NCFE 13 3 IMPLEMENTATION

Primary fields in layered media

We compute the primary fields for a media with an arbitrary number of horizontal and isotropic layers using a fortran code where the dipole point sources can be placed within any layer of the model. The primary electric and magnetic fields are calculated at every position of the domain, which is necessary to introduce them as sources of the secondary fields. When obtaining the sources for the secondary fields, it is assumed that the primary electric field is constant within each element.

The input parameters defining the 1D model are: number of layers, their conductivities and thicknesses; type, position, operating frequency, and orientation of the source; and a grid indicating the

x, y, z positions for the calculation of the field. These coordinates are the midpoints of each one of the elements in which the computational domain is divided. The computation through the grid is done for each x, y along z. This allows an efficient resolution by calculating a single x and y dependent reflection and transmission response within each layer for several values of z within the layer.

In order to validate our implementation, we first compute the solution for an homogeneous half-space where source and receiver are placed on the surface, and compare the results to the ones obtained with analytical expressions [START_REF] Ward | Electromagnetic theory for geophysical applications[END_REF]. Fig. 1 shows the results obtained for a homogeneous half-space model with conductivities σ air = 10 -14 S/m and σ b = 0.01 S/m, using the three source types (HED, HMD and VMD), and a frequency range of 10 -2 to 5 × 10 3 Hz. The horizontal dipole sources are x-directed and the VMD is z-directed. The source-receiver offset is 100 m on the x-axis.

We observed that the semi-analytical and analytical fields are practically the same for all sources in the frequency range considered. To test our code, we compared our results with those obtained using the EMDPLER algorithm [START_REF] Singh | EMDPLER: A F77 program for modelling de EM response of dipolar sources over the non-magnetic layer earth models[END_REF]. This program only allows to evaluate the fields above the surface of the 1D layered model and three types of sources are implemented (HED, HMD and VMD). In Fig. 2 we show the results for a three-layer model, with conductivities σ air = 10 -14 S/m (upper half-space), σ 1 = 0.01 S/m (500 m thick surface layer), and σ 2 = 0.1 S/m (lower half-space),

where sources and receivers are placed at z = 0 m and the source-receiver offset is 100 m on the xaxis. The horizontal and vertical sources are xand z-directed, respectively. We observe an excellent agreement between solutions in the frequency range considered.

Sparse direct solver for global finite element approach

Direct solvers have been successfully implemented for modeling (e.g. [START_REF] Streich | 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy[END_REF][START_REF] Da Silva | A finite element multifrontal method for 3D CSEM modeling in the frequency domain[END_REF][START_REF] Chung | Three-dimensional modelling of controlled-source electromagnetic surveys using an edge finite-element method with a direct solver[END_REF] and inversion (e.g. [START_REF] Grayver | Three-dimensional parallel distributed inversion of CSEM data using a direct forward solver[END_REF][START_REF] Manassero | A Reduced Order Approach for Probabilistic Inversions of 3D magnetotelluric data I: General Formulation[END_REF][START_REF] Manassero | A Reduced Order Approach for Probabilistic Inversions of 3D Magnetotelluric Data II: Joint inversion of MT and Surface-Wave Data[END_REF].

This has been possible due to the development of efficient and scalable sparse matrix factorization algorithms [START_REF] Da Silva | A finite element multifrontal method for 3D CSEM modeling in the frequency domain[END_REF]. These methods are especially efficient for multi-source problems since once the computationally demanding matrix factorization is completed, the multiple solutions are quickly computed with low computational cost [START_REF] Puzyrev | Evaluation of parallel direct sparse linear solvers in electromagnetic geophysical problems[END_REF].

Here we solve Eq. ( 7) using the package MUlti-frontal Massively Parallel Sparse direct Solver (MUMPS), version 5.1.2 [START_REF] Amestoy | MUMPS: A General Purpose Distributed Memory Sparse Solver[END_REF]. MUMPS solves the linear system on distributed memory computers by implementing a method based on a multi-frontal approximation that performs a Gaussian factorization. The system is solved in three main levels: one level of ordering and symbolic factorization, which is crucial for optimizing memory usage [START_REF] Guermouche | Impact of reordering on the Memory of a Multifrontal Solver[END_REF] and reduce the number of operations in the process of factorization and solution calculation. The next level is the factorization, which in our case is of type K = LDL T , since K is symmetric. This stage is the most demanding in terms of computational cost, but for systems with several right-hand sides only one factorization is required. The third and final stage involves solving the system through successive forward and backward elimination steps.

Parallel iterative solver for domain decomposed finite element approach

The straightforward implementation of the DDFE as presented in section 2.3 results in a slowly convergent algorithm. As demonstrated in [START_REF] Gauzellino | Nonconforming finite element methods for the three dimensional Helmholtz equation: iterative domain decomposition or global solution?[END_REF], increasing the size of the domains to include several finite elements, and solving within each of them a global problem, significantly reduces the time needed for the algorithm to reach convergence. In our implementation, the subdomains of the DD are strips in the x direction, composed of n x elements Ω j . With this implementation only n y n z linear systems of 10n x + 2 unknowns are solved instead of n x n y n z linear systems of 12 unknowns, as suggested in section 2.3. Notice that the linear system of equations in each domain does not need a parallel solver, because it is small in size; thus they are solved using a direct solver for symmetric systems. Moreover, the structure of the domain decomposition in strips improves the efficiency of the iterative process because instead of solving the linear systems and updating the Lagrange multipliers sequentially in each subdomain, it uses a blue-black iterative scheme [START_REF] Zyserman | Parallel finite element algorithm with domain decomposition for three dimensional magnetotelluric modelling[END_REF]. This scheme consists of dividing the subdomains into two types, blue and black, such that the blue ones are surrounded by the black ones (Fig. 3). The resolution of the linear system of equations and the updating of the Lagrange multipliers is carried out first for the blue, and then for the black subdomains. This enables a reduction of about 50% in the number of iterations to reach the given tolerance (Zyserman & Santos, 2000).

We use Message Passing Interface (MPI) for the parallel code implementation [START_REF] Pacheco | Parallel programming with MPI[END_REF]. We assign a portion of our domain comprised of several subdomains (strips) to each processor (which receives approximately the same number of stripes) (Fig. 3). In each iteration, each processor solves 3-D CSEM modeling using NCFE 15 the differential problem in the assigned region. The flow of information between processors sharing virtual faces of the domain ensures that local solutions converge to the global solution.

Here we evaluate the performance of nonconforming FE algorithms for modeling the 3-D CSEM response using synthetic models. We also asses their efficiency when working on distributed memory computers and their sensitivity in complex 3-D conductivity structure scenarios.

Comparison with disk models

In order to evaluate the accuracy of our DDFE and GFE algorithms, we consider two resistive-disk models representative of 3-D CSEM case studies in ocean and on land (Fig. 4), proposed by [START_REF] Constable | Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling[END_REF] and Streich & Becken (2011a), respectively.

The canonical disk model proposed by [START_REF] Constable | Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling[END_REF] (C-W disk model) represents a case study of deepwater marine CSEM and has been widely used as a testing tool in 3-D CSEM modeling (e.g. [START_REF] Streich | 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy[END_REF][START_REF] Schwarzbach | Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics-a marine CSEM example[END_REF][START_REF] Puzyrev | A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling[END_REF]. The target is a hydrocarbon reservoir 1 km below the seabed. The proposed electrical conductivity distribution consists of two halfspaces, the sea water (3.33 S/m) and underlying sediments (1.0 S/m). A 100 m thick resistive disk with conductivity σ = 0.01 S/m is placed 1 km below the interface, representing the hydrocarbon-saturated sediment region (Fig. 4a). Two different disk diameters, 2 km (small) and 5 km (large) were contemplated. In order to appreciate the influence of the airwave on the EM data [START_REF] Schwarzbach | Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics-a marine CSEM example[END_REF], we also consider a variant of this model that consists of replacing the water half-space with a 1 km layer of water and, above it, a very low conductivity air half-space (10 -14 S/m). The transmitter, an

x-directed HED, is located directly above the left edge of the resistive disk at 100 m above the seafloor and operates at a frequency of 1 Hz.

The model discretization is determined by the skin depth δ = 503 1/(σf ), where f is the operating frequency of the source and σ is the conductivity of the water or sediments which determine δ min and δ max , respectively. We use δ min to define the maximum element dimensions and 5δ max to define the minimum size of the computational domain. For the original C-W model we construct a mesh of 96 × 64 × 24 elements for the small disk case, and of 144 × 96 × 24 elements, for the large disk case.

The dimensions of the computational domain are 13 × 7.5 × 3.15 km 3 and 13 × 11 × 3.15 km 3 respectively. When considering the version that includes the air layer, the number of elements increases to 32 and the size of the domain to 4.4 km in the z-direction. As shown in Fig. 5, the mesh spacing is not uniform, with maximum values of 300 m and minimum values of 25 m, in order to achieve a proper representation of the conductivity structure and the variations of the EM field.

As usual in marine CSEM studies, the horizontal component of the electric field on the seabed is considered. In this particular case, we take from our simulated data the in-line x-component of the electric field, E x , and compute the amplitude and phase. As shown in Figs. 6 and 7 both solutions (DDFE and GFE) are practically identical. The amplitude is in excellent agreement with the synthetic data from [START_REF] Constable | Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling[END_REF], whereas the phase presents small differences. This may be due to the phase differences between the electric fields calculated in the 1D models or the different approximation strategies. The EM response to 3-D conductivity variations is satisfactory and presents no inconsistencies between the domain decomposed iterative version and global solution.

It is important to note that in order to carry out the simulation of the model containing the very low conductivity layer (air layer), we have not extended the computational domain and have kept the edges of the domain relatively close to the anomalous conductivity body. The airwave practically does not decay over the entire lateral extent of the domain and it can still be captured without response distortions at the edges of the computational domain. This can be observed in Fig. 8, where we present amplitude and phase of E x in-line models with and without the air layer. These results reproduce those presented by [START_REF] Schwarzbach | Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics-a marine CSEM example[END_REF]. The influence of the airwave becomes predominant beyond 6 to 8 km masking the anomalous response of interest.

The model proposed by Streich & Becken (2011a) (S-B disk model) is a simulation of an on land CSEM survey. The electrical conductivity model is based on a pilot CO 2 capture site in Ketzin, Germany (Streich & Becken, 2011a). The 3-D body is a disk of 2 km diameter, 15 m thickness, and 0.1 S/m conductivity, buried at 635 m depth between two sediment layers of conductivity 1/3 S/m and 1 S/m (Fig. 4b). The source is a 1 km long x-directed wire placed 0.1 m below the surface and centered directly above the left edge of the disk. The primary fields were obtained following the formulation in Streich & Becken (2011a) and applying the digital filter strategy to the expressions of the field generated by the finite length wire. The receivers are located 0.05 m below the surface.

The dimensions of the computational domain used are 7.6 × 6.6 × 2.4 km 3 and we use a mesh of 106 × 96 × 32 elements. The air layer is 600 m thick, we center the disk at (2, 3.3) km and the source at (1, 3.3) km. Comparing the electric field amplitude obtained along x-direction (in the line of the source) with the synthetic data from Streich & Becken (2011a), we observe (Fig. 9) a good agreement (both in amplitude and in location) of the anomaly. In order to show the small variations of less than 0.3% between solutions, the total field amplitude, E x (total), in Fig. 9 is normalized with respect to the primary field.

Scalability test

3-D CSEM modeling is a computationally demanding process where factors such as the size of the computational domain, the discretization method, and the numerical resolution technique become crucial to obtain solutions in practical times. In order to evaluate the scalability of our algorithms, it is necessary to test their behavior when working on distributed memory computers.

With this purpose in mind, we performed the simulations with different number of processors working in parallel: 8, 16, 32 and 64. A single node with four AMD Opteron 6276 CPUs with 16 cores each was used. We considered a reference model consisting of three layers (air half-space, 1 km thick water layer and sediments half-space with conductivities 10 -14 S/m, 3.6 S/m and 0.5 S/m, respectively)

and an anomalous conductivity layer (0.01 S/m) embedded in the sediment layer at 1 km depth from the seabed. The dimensions of the computational domain are 15 × 15 × 5 km 3 and we use different discretization sizes: 64 × 64 × 24, 96 × 96 × 24, 128 × 128 × 32 and 160 × 160 × 32 elements, corresponding to the degrees of freedom (dof) indicated in Table 2.

We evaluated the performance of our algorithms considering the total computational time without parallelization (T s ) and for N processors (T N ). We define the speed up as:

Speed up = T s T N . (14)
We compute the speed up for the different mesh sizes taking as a starting point for all: speed up = 8 for 8 processors (i.e. T s = 8T 8 ). As a reference, the ideal speed up curve is included in Fig. 10. We observe that the computational time decreases as the number of processors increases for both the iterative and direct solver. Speed up gets closer to the ideal one when the mesh size increases. On the other hand, the parallelization is more efficient for the iterative solver than for the direct solver. This results are also shown in the Table 3, where, for example, the speed up with 64 processors is approximately 30 for DDFE while the maximum value is 20.1 for the GFE. In addition, we can appreciate that the computation times between algorithms are considerably different and do not increase linearly as the mesh size increases. This is important when designing realistic models according to our requirements and computational resources.

The computation time of the algorithms can be dependent on the model's conductivity values and the source's operating frequency. Based on a benchmark exercise by [START_REF] Weiss | Project APhiD: A Lorenz-Gauged A-Φ Decomposition for Parallelized Computation of Ultra-Broadband Electromagnetic Induction in a Fully Heterogeneous Earth[END_REF], we evaluate our algorithms' performances according to the aforementioned parameters. Two models are considered, which are slightly different from the model described above. The first one has a water half-space ("deep water" model), and the second one has an air half-space and a 500 m thick water layer ("shallow water" model). We also consider frequencies of 10, 1, and 0.1 Hz, and the grids of 64 × 64 × 24 and 128 × 128 × 32 elements to make comparisons. The DDFE algorithm presents a monotonically decreasing relative error, smooth convergence, and rapidly reaches the set tolerance of 10 -4 for the relative error in a few hundreds iterations (Fig. 11). We note that there are no significant changes in computation time, nor number of iterations, due to the presence of the very low conductivity layer in the "shallow water" model. However, this is not the case with respect to frequency, where faster convergence is observed for the higher frequencies. In Table 4 the computation times are compared 3-D CSEM modeling using NCFE 19

for the GFE algorithm, where it is shown that no significant differences were observed with respect to models and frequency.

Real-world conductivity structure model

Here we analyze the performance of the algorithms when dealing with complex geological structures.

At the same time we want to evaluate the sensitivity of the method to anomalous conductivity bodies in a realistic and challenging setup. For this, we propose a 3-D synthetic model based on a real CSEM case study offshore Brazil developed by YPF S A.

The model dimensions are 30 × 18 × 9.9 km 3 . The different layers (air, sea water and sediments) and their respective conductivities are shown in Fig. 12. The geometry is based on the seismic acquisition in the area, and the conductivity values are based on the geological interpretations of the area. We consider a flat seafloor located at z = -2.7 km. The source is a HED working at 1, 0.5, 0.25 and 0.05 Hz at 50 m above the seafloor. The source is set in six different positions along two lines, L1 and L2 (Figs 12d and e, and Table 5). In each case, the source is oriented in the line direction. A mesh of 96 × 64 × 56 elements with varying elements sizes was used to obtain an accurate representation of the structures throughout the domain.

We place the oil reservoir, to which we assign a bulk electrical conductivity of 0.05 S/m, at different locations within the second layer of sediments, which in turn has a conductivity of 1.6 S/m. The proposed sizes and locations of the reservoir are shown in Figs 12(a), (b) and (c), and referred to as case 1, 2 and 3, respectively. We assume that the depth of the layers have been previously determined with data from seismic reflection and wells, and use this background model as a reference for the EM response. In order to separate the response of the reservoir, we normalize the full EM response with the response of the background model. Note that this representation does not eliminate the dependence of the field amplitude with the source position.

We use both algorithms (DDFE and GFE) in series to evaluate the CSEM response at the four chosen frequencies and the six different source positions. The survey design implies accounting for different EM sources at different positions and directions. In this case, the advantage of obtaining direct solutions is remarkable, since the analysis and factorization of the matrix of Eq. 7 (∼ 98% of the computational time) is done only once and is reused for each source. Thus, the total computational time (solution for 4 frequencies and 6 sources) taken by the DDFE algorithm was 27.0 min, while the GFE algorithm runs in 20.6 min.

As usual in marine CSEM surveys, we are interested in the horizontal electric field at the seafloor level to interpret our data. Since the marine CSEM method is strongly dependent on the source-receiver configuration with respect to the reservoir [START_REF] Um | On the physics of the marine controlled-source electromagnetic method[END_REF], we expect to determine preferential acquisitions where it would be most evident to locate the anomaly corresponding to the presence of an oil reservoir. We find that the highest amplitude of the field anomalies are obtained when the source is laterally shifted with respect to the reservoir's location. In Fig. 13 we present the amplitude of E x normalized with respect to the reference field for a frequency of 0.5 Hz. In all cases, the responses of the model with the reservoir are 10% higher than the reference field, and there is a spatial variation of the anomaly associated to the lateral conductivity distribution. Since the field is normalized, the amplitude of the anomaly is affected by the amplitude of the reference field. This can be observed in case 3 with source 3 (Fig. 13 left bottom), where the normalized field increases as it approaches an angle of approximately 45 degrees to the x-directed HED. Along these directions the amplitude of the

x-component of the primary field is minimum. In Fig. 14 we show the normalized amplitude of E x in-line for the multiple frequencies considered. On line 1, in all three cases, it can be seen that the field is sensitive to lateral variations in conductivity, especially on the edge of the reservoir located farthest away from the source. Comparing cases 2 and 3, with source 3 or 6, we can see significant changes in amplitude, even if there is only a 500 m difference in the depths of each reservoir. We note that, in this particular case, it seems preferable to use high frequency, since this produces greater amplitude than lower frequencies.

We have introduced and tested the performance of two algorithms (DDFE and GFE), which solve the CSEM response of 3-D media using nonconforming finite element methods. Both algorithms were shown to accurately reproduce the CSEM response of canonical models, and show competent speed up curves. In addition, for all the models presented in this work, the numerical codes can be employed on personal computers. This is a valuable quality since large clusters are not always accessible to the scientific community.

Although we have not found definite advantages of one algorithm over the other, they do have differences that make them more suitable than the other in specific applications. The code developed with the DDFE algorithm is self-contained, that is, it does not require installation of external libraries, and its use is immediate. In contrast, the GFE implementation requires the MUMPS package and multiple associated libraries. Therefore, its installation can be cumbersome and different compatibility problems can appear. This could be an inconvenient for researchers working in groups without system administrator's support. In particular, we observed that the efficiency of MUMPS, and therefore of the GFE implementation, is not the best if the BLAS libraries available are not the optimized ones. To get them, some skill beyond scientific programming is needed. In our real-world conductivity structure example, we have seen that the use of the GFE algorithm might be advantageous when dealing with multi-source problems because it efficiently solves systems with multiple right-hand sides.

The structured grid used in our approach could be considered in principle not flexible enough. But given the presented speed up curves and possibilities of optimization in terms of multi-source simulations, we find that these grids are a suitable alternative given the simplicity of their implementation.

Moreover, the good results we have obtained in terms of computing times naturally lead us to think of incorporating our numerical forward codes into probabilistic inversion algorithms to estimate the subsurface conductivity from CSEM data. In such strategies, it is standard practice to divide the 3-D domain of investigation in regular voxels and use simple neighbor-based regularization schemes to constrain the properties of the unknown variables [START_REF] Rosas-Carbajal | Two-dimensional probabilistic inversion of plane-wave electromagnetic d ata: Methodology, model constraints and joint inversion with electrical resist ivity data[END_REF][START_REF] Rosas-Carbajal | Probabilistic three-dimensional time-lapse inversion of magnetotelluric data to infer mass transfer in a geothermal system[END_REF][START_REF] Manassero | A Reduced Order Approach for Probabilistic Inversions of 3D magnetotelluric data I: General Formulation[END_REF].

Our 3-D numerical codes allowed us to evaluate the potential contributions of acquiring CSEM data in a realistic offshore scenario. The 3-D geometry and different options for the reservoir locations (Fig. 12) were built based on real seismic data and their interpreted sections. Our results show that the presence of an oil reservoir in any of the three cases proposed would create a response in the data significantly different from the response without a reservoir. Note that this exercise was done assuming a perfect knowledge of the subsurface geometry, and does not mean that the exact location of the anomaly would be easy to determine. However, it is highly encouraging because it shows that a CSEM survey could help to determine which of the 3 reservoir scenarios is most likely. More specifically, the visualization of the anomalous response along line 1 (Fig. 14) allows us to optimize the relationships between frequency, amplitude and target depth. We see that whenever a symmetry axis is identified, it is easier to recognize and define the location of the anomaly with a survey line oriented perpendicular to these symmetry axis.

The formerly discussed capabilities of our algorithms naturally carry our future work into the design of an efficient data inversion scheme. For this, it will be necessary to evaluate the performance of the algorithms in their frequency parallelization, that is, determining which algorithm is more convenient when distributing to different processors the tasks of solving the same problem for different frequencies. Due to the way the solution is computed, we expect the GFE algorithm to be better disposed to this type of parallelization, since the factorization time depends mainly on the problem size, and not on the frequency. We expect that this strategy, coupled with the efficient solution of multi-source problems, will be the most efficient way to implement a 3-D CSEM inversion algorithm.

We have presented two 3-D modeling algorithms using nonconforming finite elements with 1) a direct solver (GFE), and 2) an iterative solver (DDFE), to simulate the CSEM response in the frequency domain operating with electric or magnetic dipole sources. Simulations were carried out to evaluate the computational performance of the algorithms and the accuracy of the solutions.

The primary electromagnetic field generated by dipole electromagnetic sources in a layered media is particularly useful in CSEM to accurately model the field in the source vicinity and the presence of the air layer in the models. The code computing the primary field solution has been compared with analytical solutions and those yielded by the EMDPLER code solution with excellent agreement over a wide frequency range. Also, the comparison of our solutions obtained using the NCFE algorithms with those of other approaches described in the literature are in excellent agreement.

We used a scalability analysis to evaluate the parallelization performance of the algorithms. We found that the parallelization is more efficient when the mesh size increases, and that the DDFE is better suited to parallelization than the GFE. In addition, our codes offer satisfactory parallel efficiency with a considerable reduction in total computation time as the number of processors increases.

Finally, we have shown the capability of our 3-D CSEM modeling scheme to deal with realistic geological scenarios considering multiple sources and frequencies. The modeling was done inspired in a real offshore seismic and CSEM acquisition, and we showed that even in this complex 3-D geometry the EM data could help to determine whether an oil reservoir is located within a sedimentary layer. with synthetic data from [START_REF] Constable | Marine electromagnetic methods-A new tool for offshore exploration[END_REF]. DDFE and GFE solutions fit the reference data although with small differences for source-receiver offsets greater than 1 km. 3-D CSEM modeling using NCFE 37

Figure 1 .

 1 Figure 1. Comparison between the computed primary field solution and those obtained from analytical expressions. Three types of dipole sources are presented: HED, HMD and VMD. The source and receiver are located at z = 0 m in a homogeneous half-space of conductivities σ air = 10 -14 S/m and σ b = 0.01 S/m. The sourcereceiver offset is 100 m. The curves are in agreement over the frequency range usually used in CSEM.

Figure 2 .

 2 Figure 2. Comparison between the computed primary field solution and those obtained with the EMDPLER code. The three types of sources implemented in the EMDPLER code are presented. We simulate a three-layer model of conductivities σ air = 10 -14 S/m (upper half-space), σ 1 = 0.01 S/m (500 m thick surface layer), and σ 2 = 0.1 S/m (lower half-space). The source and receiver are located at the interface between the air layer and the first layer. The source-receiver offset is 100 m.

Figure 3 .

 3 Figure 3. Left: the blue-black scheme used to speed up the iterative process. Each strip (domain) involves n x elements of the finite element grid. Right: scheme of the division of the domain Ω assigned to each processor. The arrows indicate the flow of information between adjacent cells.

Figure 4 .

 4 Figure 4. Vertical slice of the resistive-disk models (a) offshore, proposed by Constable & Weiss (2006), and (b) onshore, proposed by Streich & Becken (2011a). In each case, the vertical plane exactly slices the disk in half along the x-direction and the dimensions of the computational domain and the disk are only schematic.

Figure 5 .

 5 Figure 5. (a) and (b) Slices of the mesh used for the discretization of the Constable & Weiss (2006) model with the 2 km diameter disk, (a) in the plane y = 3.75 km and (b) in the plane xy at z = 2 km. (c) Slice of the mesh used for the version of the C-W model, that includes the air layer, in the plane y = 3.75 km.

Figure 6 .

 6 Figure 6. Comparison between our solutions and the synthetic data from Constable & Weiss (2006). The amplitude of the electric field in-line x-component, at z = 1 km depth (at the seafloor). Both solutions, DDFE and GFE, are identical to each other and are in excellent agreement with the synthetic data presented by Constable & Weiss (2006).

Figure 7 .

 7 Figure 7. Comparison of the phase of the electric field in-line x-component, at z = 1 km depth (at the seafloor),

Figure 8 .

 8 Figure 8. Amplitude and phase of the electric field in-line x-component of the Constable & Weiss (2006) model with and without air layer, at z = 1 km (at the seafloor). The influence of the air layer is seen to dominate the EM response at source-receiver offset greater than 6 km.

Figure 9 .

 9 Figure 9. Comparison between the DDFE and GFE solutions and the synthetic data of Streich & Becken (2011a). Normalized in-line x-component electric field with respect to the background field, at z = 0.05 m below the surface. Small differences in amplitude less than 0.3% of the background field are observed.

Figure 10 .

 10 Figure 10. Speed up curves of the DDFE and GFE algorithms using different grid sizes. The computation was carried out with 8, 16, 32 and 64 processors. The speed up achieved remains relatively constant.

Figure 11 .

 11 Figure 11. Convergence curves of the DDFE algorithm for two models, "deep water" and "shallow water", using different frequencies and grid sizes. The presence of the air represented by the very low conductivity layer in the "shallow water" model does not affect convergence time, whereas, clearly, the frequency and grid size does.

Figure 12 .

 12 Figure 12. Realistic synthetic model. (a), (b), and (c) vertical slices of the different reservoir positions, case 1, 2, and 3, respectively. (d) Horizontal slice at z = -5 km of case 1. Lines L1 and L2 represent the transects, S1 to S6 denote the position of the different sources. (e) 3-D view of the conductivity model of case 1 with the plotting of lines L1 and L2.

Figure 13 .

 13 Figure 13. x-component of the electric field (E x (with reservoir)) normalized with respect to the background field (E x (without reservoir)), at frequency f = 0.5 Hz. The fields are measured at the seafloor, located at z = -2.7 km. Top, medium and bottom rows correspond to cases 1, 2 and 3, respectively (see Fig. 12 a, b and c). Left column shows sources along L1, right one along L2. In each case, the anomalous response is located around the reservoir location, reaching amplitudes greater than 10% of the background field.

Figure 14 .

 14 Figure 14. x-component of the electric field (E x (with)) normalized with respect to the background field (E x (without)), at frequencies f = 1, 0.5, 0.25, and 0.05 Hz. The fields are measured at the seafloor, located at z = -2.7 km. Top, medium and bottom rows correspond to cases 1, 2 and 3, respectively (see Fig. 12 a, b and c). Left column shows sources along L1, right one along L2. In each case, the relation between frequency and amplitude of the EM response and depth of the target can be observed.

Table 4 .

 4 GFE algorithm computation times for two models, "deep water" and "shallow walter", using different frequencies and grid dimensions.

			Deep water		Shallow water
	Frequency [Hz]	10.0	1.0	0.1	10.0	1.0	0.1
	Grid dimension 64 × 64 × 24
	CPU time [s]	40.5	41.8	42.1	48.3	51.5	51.0
	Grid dimension 128 × 128 × 32
	CPU time [s]	458.9 458.3 453.7	379.0 379.0 485.7

Table 5 .

 5 Source position on each line. The HED is located 50 m above the seabed and oriented in the direction of the corresponding line.

	Line 1 (x,y) [km] Line 2 (x,y) [km]
	Source 1 (10,9)	Source 4 (8,16)
	Source 2 (15,9)	Source 5 (15,9)
	Source 3 (20,9)	Source 6 (22,2)

ACKNOWLEDGMENTS

We wish to thank Dr. R. Streich for providing the data used in the comparison of Fig. 9. We thank Geoph. Federico G. E. Späth (YPF SA) for the exchange of ideas and important discussions for the real-world model in Section 4.3.

The numerical simulations were performed in TUPAC, the CSC-CONICET computing cluster; special thanks to the system administrators D. Vinazza and E. Mocksos for their helpfulness. Finally, we thank the constructive comments by reviewers J. D. De Basabe and C. Weiss, which contributed to improve this work.

DATA AVAILABILITY

The source code and all data produced in this work are available and can be requested to the e-mail melias@fcaglp.unlp.edu.ar.

Table 1. Basis functions used for the nonconforming finite element method. Modified from Zyserman & Santos (2000).