# Anisotropic three-dimensional inversion of CSEM data using finite-element techniques on unstructured grids

Anisotropic three-dimensional inversion of CSEM data using finite-element techniques on... Summary In this paper, we present a recently developed anisotropic 3-D inversion framework for interpreting controlled-source electromagnetic (CSEM) data in the frequency domain. The framework integrates a high-order finite-element forward operator and a Gauss–Newton inversion algorithm. Conductivity constraints are applied using a parameter transformation. We discretize the continuous forward and inverse problems on unstructured grids for a flexible treatment of arbitrarily complex geometries. Moreover, an unstructured mesh is more desirable in comparison to a single rectilinear mesh for multisource problems because local grid refinement will not significantly influence the mesh density outside the region of interest. The non-uniform spatial discretization facilitates parametrization of the inversion domain at a suitable scale. For a rapid simulation of multisource EM data, we opt to use a parallel direct solver. We further accelerate the inversion process by decomposing the entire data set into subsets with respect to frequencies (and transmitters if memory requirement is affordable). The computational tasks associated with each data subset are distributed to different processes and run in parallel. We validate the scheme using a synthetic marine CSEM model with rough bathymetry, and finally, apply it to an industrial-size 3-D data set from the Troll field oil province in the North Sea acquired in 2008 to examine its robustness and practical applicability. Electrical anisotropy, Marine electromagnetics, Inverse theory, Numerical approximations and analysis 1 INTRODUCTION For derisking the exploration of submarine oil and gas resources, the controlled-source electromagnetic (CSEM) method is now regarded as a crucial surveying methodology providing electrical resistivity information on the subsurface (Kong et al.2002; Johansen et al.2005; Constable & Srnka 2007; Holten et al.2009; Anderson & Mattson 2010; Schwalenberg & Engels 2012; Fanavoll et al.2014a,b; Hölz et al.2015). It is complementary to seismics and emerged in the 2000s (Eidesmo et al.2002). A horizontal electric dipole (HED) source can be deployed efficiently and generates a field with feasible polarization to detect thin resistive hydrocarbon-bearing layers. Nowadays, the advances in the development of instruments have greatly improved the data quality and acquisition efficiency such that 3-D dense survey layouts are becoming routine for acquiring large volumes of data (Gabrielsen et al.2009; Pedersen & Hiner 2014). A comprehensive data coverage can enhance the imaging resolution of 3-D structures and thus reduce the ambiguity in data analysis. However, this also creates the need for developing sophisticated interpretation tools such as 3-D inversion algorithms to treat large data sets. An important aspect of interpreting marine CSEM data is to understand the effect of electrical anisotropy (Løseth et al.2007; Ellis et al.2010). The fine layering of marine sediments and asymmetric grain shapes yield a larger resistivity in the direction normal to the deposited layers than parallel to them (Tompkins 2005; Newman et al.2010; Brown et al.2012). In this paper, we assume that the layers have small tilt angles so that the layer normal and parallel directions can be approximated as the vertical and horizontal directions (VTI anisotropic medium, Jing et al.2008). The impact of electrical anisotropy on CSEM data has a complicated dependence on many factors such as water depths and survey geometries. Based on a 1-D scenario with varying degrees of electrical anisotropy, Lu & Xia (2007) illustrated that the inline data are stronger related to the vertical resistivity, while the broadside measurements are more sensitive to the horizontal resistivity. This observation was further elaborated by Ramananjaona et al. (2011) using 1-D forward modeling and inversion. Furthermore, they also observed that inversion of marine CSEM data measured in an anisotropic environment under the assumption of electrical isotropy can lead to a resolved resistivity model with many misleading oscillations (artefacts). Similar conclusions have also been found in both 2-D (Ramananjaona & MacGregor 2010) and 3-D inversion scenarios (Jing et al.2008; Newman et al.2010; Sasaki 2011). Although broadside data are less sensitive to thin resistive layers than inline data, the incorporation of them in an anisotropic inversion is proven to better resolve the electrical anisotropy (Newman et al.2010) and to be useful to reduce the overall data misfit (Mohamad et al.2010). These findings have been promoting the use of full 3-D data sets. Besides allowing for the host medium being electrically anisotropic, Abubakar et al. (2010) and Brown et al. (2012) have discussed the possibility of resolving the electrical anisotropy of a thin reservoir, proposing the addition of other sources such as horizontal magnetic dipoles to enhance the inductive effect for reconstructing the horizontal resistivity. In this paper, however, we only focus on an HED source due to the field data set at hand. In order to interpret marine CSEM data, reported anisotropic 3-D inversion schemes include a parallel non-linear conjugate gradient (CG) optimization approach with a finite-difference forward operator (Newman et al.2010), a quasi-Newton optimizer equipped with a time-domain finite-difference forward code (Zach et al.2008), a Gauss–Newton optimization with a finite-difference forward operator (Sasaki 2013; Amaya et al.2016), a reweighted regularized CG method with an integral equation modeling operator (Zhdanov et al.2014b) and a Schur complement finite-difference forward operator (Jaysaval et al.2014) for serving the anisotropic quasi-Newton optimization approach presented by Zach et al. (2008). In this paper, we introduce a new anisotropic 3-D inversion code which has been implemented using finite-element techniques. Although finite-element methods have received considerable attention for electromagnetic forward modeling (Schwarzbach et al.2011; Ren et al.2013; Ansari & Farquharson 2014; Grayver & Bürg 2014), to our best knowledge, the published inversion results utilizing finite-element approaches were mainly 2-D anisotropic (Key 2012; Key et al.2014) and 3-D isotropic (Schwarzbach & Haber 2013; Grayver 2015; Zhang & Key 2016). Rather than discretizing the forward problem and parametrizing the inverse problem on structured grids, we adopt unstructured mesh discretization. Unstructured meshes offer great flexibility in handling complex geometries to allow for very accurate descriptions of complicated seafloors, reservoir structures, or submetre scale infrastructure like pipelines and casings in the geomodel. Unstructured meshes also make it extremely convenient to carry out local grid refinement. High-order finite-element methods are available for solving the forward problem with high accuracy (Castillo et al.2005; Schwarzbach 2009). We formulate the forward problem by means of the so-called secondary field approach which only simulates scattered fields generated by the discrepancies of electromagnetic parameters from a known background for primary fields (Streich 2009; Schwarzbach et al.2011). Although the total field approach is also implemented in the code, the secondary field approach is favoured here for simulating marine CSEM data. The reason for this is twofold. At first, it is reasonable to assume that the sea surface is flat such that a horizontally stratified model is readily regarded as the background in which field responses can be efficiently obtained. Moreover, the singularity removal technique that the secondary field method applies allows a reduction of the node density around the sources which makes it advantageous over the total field approach when many sources of finite length are considered in the same model (mesh). As marine CSEM data are often processed at only a few frequencies but many transmitter positions (Gabrielsen et al.2009; Weitemeyer et al.2010; Zhdanov et al.2014a), it is desirable to use a parallel direct solver to inexpensively obtain the responses of the transmitters operated at the same frequency. We formulate the inverse problem in an unconstrained least-squares sense, penalized by the smoothness regularization. Although it has been argued that images stabilized by smoothness functionals cannot sufficiently describe the true geological structures (Portniaguine & Zhdanov 1999; Gribenko & Zhdanov 2007; Zhdanov et al.2007), our experiences with this regularization are positive because it is very robust and barely problem-dependent (such as on model parametrization). We solve the discrete inverse problem using an inexact Gauss–Newton approach by virtue of its relatively fast convergence rate. The scheme is implemented in a parallel manner to accelerate the inversion of multisource and multifrequency EM data. Its efficiency and versatility are further considered by providing three alternatives to handle the Jacobian. In the following, we first revisit the finite-element discretization of the forward problem in terms of the secondary electric field. Then, we describe the inverse problem and introduce the inexact Gauss–Newton algorithm. To compute the Jacobian occurring in the Gauss–Newton equation systems, we offer three approaches. The efficiency of each method depends on the survey geometry and problem size and, therefore, the selection is data-driven. Finally, we examine the feasibility of the framework using a synthetic marine example with rough bathymetry and the Troll 3-D CSEM field data set acquired in 2008 in the North Sea (Gabrielsen et al.2009). 2 FORWARD PROBLEM We assume a harmonic time dependence e−iωt to simulate EM fields in the space $$\mathbf {r}\in \mathbb {R}^3$$ induced by an electric source with a current density Js = Js(r, ω). We construct a bounded computational domain $$\Omega \in \mathbb {R}^3$$ with its outer boundary $$\partial \Omega \in \mathbb {R}^3$$ far enough from the source to properly characterize the true field behaviour. The associated boundary value problem consists of the curl–curl equation of the electric field and a simple Dirichlet boundary condition imposed at the outer boundary ∂Ω to give a unique solution   \begin{eqnarray} \nabla \times \,({{{\mu }}}^{-1}\nabla \times \,{\mathbf {E}}) - i \omega \boldsymbol{\tilde{\varsigma }}{\mathbf {E}}= i \omega \mathbf {J}_s \quad {\rm in} \quad \Omega , \end{eqnarray} (1a)  \begin{eqnarray} \mathbf {n}\times {\mathbf {E}}= \mathbf {0} \quad {\rm on} \quad \partial \Omega , \end{eqnarray} (1b) where n denotes the unit normal vector to the outer boundary, μ is the magnetic permeability which is assumed to be isotropic and allowed to vary throughout the model here and $$\boldsymbol{\tilde{\varsigma }}$$ denotes the tensor complex electrical conductivity. In a medium with transverse isotropy about a vertical axis, $$\boldsymbol{\tilde{\varsigma }}$$ is characterized by a 3 × 3 diagonal matrix   $$\boldsymbol{\tilde{\varsigma }}= \left(\begin{array}{ccc}\sigma _h - i \omega \epsilon _h &0 &0 \\ 0 & \sigma _h - i \omega \epsilon _h &0 \\ 0 & 0 &\sigma _v - i \omega \epsilon _v \end{array} \right),$$ (2)where σh and σv mean the horizontal and vertical conductivities, respectively, and εh and εv denote the horizontal and vertical electrical permittivities, respectively. With the working frequencies in the range between about 0.1 and 10 Hz (Constable 2010) and by assigning a small conductivity value such as 10−7 Sm−1 to the air layer (Jaysaval et al.2014), the marine CSEM method works in the quasi-static limit where σ ≫ ωε holds. However, we preserve the full-wave expression to generalize the applicability of our code. To avoid the numerical difficulty in directly simulating the rapidly decaying field around the transmitter, the total electric field E = Ep + Es is split into a known primary field Ep due to the impressed current Js in a background model ($${{\mu }}_0, \boldsymbol{\tilde{\varsigma }}_0$$) and a secondary electric field Es to be solved for. In doing so, the mesh node density around the source can be relaxed compared to the one for formulation (1). Moving the terms including Ep to the right-hand side of eq. (1) attains the boundary value problem in terms of Es (Schwarzbach et al.2011)   \begin{eqnarray} \nabla \!\times ({{{\mu }}}^{-1}\nabla \!\times {\mathbf {E}^s}) \!-\! i \omega \boldsymbol{\tilde{\varsigma }}{\mathbf {E}^s}\!=\! \nabla \!\times (\delta {{{\mu }}}^{-1}\nabla \!\times {\mathbf {E}^p}) \!-\! i \omega \delta \boldsymbol{\tilde{\varsigma }}{\mathbf {E}^p}\ \ {\rm in } \quad \Omega ,\nonumber\\ \end{eqnarray} (3a)  \begin{eqnarray} \mathbf {n} \times {\mathbf {E}^s}= \mathbf {0} \ \ {\rm on }\quad \partial \Omega , \end{eqnarray} (3b) where the right-hand side of eq. (3a) now acts as the secondary sources and $$\delta {{{\mu }}}^{-1}= {{{\mu }}}^{-1}_0 - {{{\mu }}}^{-1}$$ and $$\delta \boldsymbol{\tilde{\varsigma }}= \boldsymbol{\tilde{\varsigma }}_0 - \boldsymbol{\tilde{\varsigma }}$$ are the parameter discrepancies between the background and total models. We solve boundary value problem (3) using the Nédélec curl-conforming finite-element method on unstructured grids (Castillo et al.2005; Schwarzbach et al.2011). The computational domain Ω is decomposed into a set of non-overlapping tetrahedra first. The unknown function Es is then approximated in the curl-conforming finite-element function space which is spanned by small-support basis functions   $$\mathbf {E}^s(\mathbf {r}) \approx \sum _{i=1}^{N_{\mathrm{dof}}} E^s_i \mathbf {N}_i(\mathbf {r}),$$ (4)where Ndof denotes the total number of degrees of freedom, Ni means the ith Nédélec curl-conforming basis function (Monk 2003), and $$E^s_i$$ denotes the ith discrete unknown coefficient in the mesh. After assembling all the finite-element matrices, solving the original boundary value problem (3) is reduced to solving the following system of linear equations   $$\mathbf {K} \mathbf {x} = \mathbf {F},$$ (5)where $$\mathbf {x} = (E_1^s, \ldots , E^s_{N_{\mathrm{dof}}})^T \in \mathbb {C}^{N_{\mathrm{dof}}}$$ denotes the solution vector, $$\mathbf {K} \in \mathbb {C}^{N_{\mathrm{dof}}\times N_{\mathrm{dof}}}$$ is the system matrix which is very sparse and symmetric, and $$\mathbf {F} \in \mathbb {C}^{N_{\mathrm{dof}}}$$ represents the secondary source vector which possesses non-zero entries where the electrical or magnetic parameters of the model differ from the background. In detail, the elements of K and F read (Schwarzbach et al.2011)   \begin{eqnarray} {K}_{ij} \!=\! \int _{\Omega } \nabla \times \,\mathbf {N}_i \cdot ({{{\mu }}}^{-1}\nabla \times \,\mathbf {N}_j) {d\Omega }- i \omega \int _{\Omega } \mathbf {N}_i \cdot \left(\boldsymbol{\tilde{\varsigma }}\mathbf {N}_j \right) {d\Omega }, \end{eqnarray} (6a)and   \begin{eqnarray} F_{i} \!=\! \int _{\Omega } \nabla \times \,\mathbf {N}_i \cdot (\delta {{{\mu }}}^{-1}\nabla \times {\mathbf {E}^p}) {d\Omega }\!-\! i \omega \int _{\Omega } \mathbf {N}_i \cdot \left(\delta \boldsymbol{\tilde{\varsigma }}{\mathbf {E}^p}\right) {d\Omega }. \end{eqnarray} (6b) Due to the sparsity of the matrix K, the system of eq. (5) can be solved either iteratively or directly. Considering that realistic marine surveys often yield thousands of transmitter positions, we prefer a parallel direct solver which might outperform an iterative approach for moderate-size problems because the solutions for the transmitters at one frequency can be computed rapidly by merely forward and backward substitutions once the system matrix is decomposed and the factors are stored. 3 INVERSE PROBLEM We discretize the forward and inverse problems on the same unstructured mesh which allows for incorporating all transmitters, frequencies, and receivers, and, at the same time, provides a reasonable inversion resolution. Assuming that the inversion domain Ωm coincides with or is a part of the forward model Ω: Ωm ⊆ Ω and includes Nm non-overlapping subdomains, we parametrize the conductivity model using the pixel-based method (e.g. Commer & Newman 2009). In each subdomain, the horizontal (σh) and vertical (σv) conductivities are decoupled and treated as independent unknown parameters   $$\left( \begin{array}{c}\sigma _h(\mathbf {r})\\ \sigma _v(\mathbf {r}) \end{array} \right) \approx \sum _{i=1}^{N_m} \psi _i(\mathbf {r}) \left( \begin{array}{c}\sigma _{h_i}\\ \sigma _{v_i} \end{array}\right),$$ (7)where $$\sigma _{h_i}$$ and $$\sigma _{v_i}$$ refer to the unknown horizontal and vertical conductivity coefficients (degrees of freedom) with respect to the ith inversion subdomain $$\Omega _{m_i}$$, respectively, and ψi(r) denotes the corresponding basis function which is non-trivial and equals unity only over $$\Omega _{m_i}\!$$. To ensure that the reconstructed electrical conductivity is always positive and finite, we enforce the bounded model transformation with a specified lower and upper conductivity bound set to σa and σb, respectively (Commer & Newman 2008). Instead of directly recovering the electrical conductivity, the transformed parameter m will be solved for first. The true conductivity is then obtained by the inverse model transformation. The forward and inverse transformations are given by   \begin{eqnarray} \left( \begin{array}{c}m_{h_i}\\ m_{v_i} \end{array} \right) \!=\!\left( \begin{array}{c}\log (\sigma _{h_i} \!-\! \sigma _a)\!-\! \log (\sigma _b \!-\! \sigma _{h_i})\\ \log (\sigma _{v_i} \!-\! \sigma _a)\!-\! \log (\sigma _b \!-\! \sigma _{v_i}) \end{array} \right); \sigma _a<\sigma _{h_i}, \sigma _{v_i} < \sigma _b,\nonumber\\ \end{eqnarray} (8a)  \begin{eqnarray} \left( \begin{array}{c}\sigma _{h_i}\\ \sigma _{v_i} \end{array}\right) = \left( \begin{array}{c}\frac{\sigma _a+\sigma _be^{m_{h_i}}}{1+e^{m_{h_i}}}\\ \frac{\sigma _a+\sigma _be^{m_{v_i}}}{1+e^{m_{v_i}}} \end{array}\right); -\infty <m_{h_i}, m_{v_i} < \infty . \end{eqnarray} (8b) It should be mentioned that the horizontal and vertical conductivity bounds (namely σa and σb in the above equations) do not necessarily have to be identical. As most typically in the VTI assumption, the horizontal conductivity is higher than the vertical one, a practical constraint for the horizontal conductivity may have a strict lower bound and a relaxed upper bound. In contrast, a loose lower bound and a tight upper bound can be more effective for the vertical conductivity. However, in the absence of any a priori information, we have used the same bounds for all the examples in the paper. Given the parametrized model and assuming that the measurements $$\mathbf {d}^\mathrm{obs}\in \mathbb {C}^{N_d}$$ are collected at Nf frequencies, Ns transmitters (at $${\mathbf {r}^s_{i}}, i = 1,\ldots ,N_s$$), and Nr receivers (at $${\mathbf {r}^r_{i}}, i = 1,\ldots , N_r$$) for Nt field components (labeled in a vector $$\boldsymbol{\Theta }$$), namely Nd = Nf × Ns × Nr × Nt, we seek a parameter model m which minimizes the following regularized discrete objective functional   \begin{eqnarray} \Phi (\mathbf {m}) &=& \Phi _d(\mathbf {m}) + \lambda _h \Phi _{m_h}(\mathbf {m})+\lambda _v \Phi _{m_v}(\mathbf {m}) \nonumber\\ &=& \frac{1}{2} \Vert \mathbf {W}(\mathbf {d}^\mathrm{obs}- \mathbf {d}^{\mathrm{syn}}(\mathbf {m}^{}))\Vert ^2 + \frac{\lambda _h}{2}\Vert \mathbf {L}_h(\mathbf {m}^{}_h-\mathbf {m}^{\mathrm{ref}}_h)\Vert ^2 \nonumber\\ &&+ \frac{\lambda _v}{2}\Vert \mathbf {L}_v(\mathbf {m}^{}_v-\mathbf {m}^{\mathrm{ref}}_v)\Vert ^2, \end{eqnarray} (9)where, in the data misfit functional, that is, the first term on the right-hand side, dsyn(m) represents the predicted data vector and W is a diagonal matrix storing the data weights. The determination of proper data weights is crucial for a successful inversion of CSEM data because the long-offset fields are orders of magnitude weaker than the short-offset ones. To treat the short- and long-offset data somehow equally, a commonly used data weight involves the standard deviation of a datum (denoted by δ|dobs| for the amplitude and δϕobs for the phase throughout the rest of the paper) (e.g. Grayver et al.2013; Schwarzbach & Haber 2013). We construct the data weighting matrix for the synthetic data set according to Schwarzbach & Haber (2013) which incorporates the consequences of the number of data, the relative error in the data amplitude, a scaling factor for different data components and the noise floor. As for the field data, a different weighting is employed to account for the rotation uncertainties which will be illustrated in Section 4.2. In the VTI inversion, the model vector $$\mathbf {m} = \left({\mathbf {m}^{}_h}^T, {\mathbf {m}^{}_v}^T\right)^T$$ is composed of the horizontal model parameter vector $$\mathbf {m}^{}_h = ({m}^{}_{h_{1}}, \ldots , {m}^{}_{h_{N_m}})^T \in \mathbb {R}^{N_m}$$ and the vertical model parameter vector $$\mathbf {m}^{}_v = ({m}^{}_{v_{1}}, \ldots , {m}^{}_{v_{N_m}})^T \in \mathbb {R}^{N_m}$$ which results in twice as many model parameters as the isotropic inversion; λh and λv in eq. (9) are the trade-off or regularization parameters; $$\Phi _{m_h}(\mathbf {m})$$ and $$\Phi _{m_v}(\mathbf {m})$$ are the model misfit functionals for the horizontal and vertical parameter models, respectively, which are chosen to be the discrete form of the smoothness stabilizer measuring the gradient of the model discrepancy in the l2 norm   $$\Phi _m(m) = \frac{1}{2}\int _{\Omega _m} \vert \nabla (m-m^{\mathrm{ref}}) \vert ^2 {d\Omega }\qquad {\rm in} \ {\Omega _m}.$$ (10)On an unstructured mesh, the discrete form of the above functional can be obtained by resorting to the lowest order divergence-conforming Nédélec finite-element method (Schwarzbach & Haber 2013). Lh and Lv in eq. (9) are identical representing the resulting smoothness matrix which is independent of the model parameters; mrefh and mrefv are the reference models for $$\mathbf {m}^{}_h$$ and $$\mathbf {m}^{}_v$$, respectively. To handle the data topology, a global indexing matrix $$\boldsymbol{\mathcal {I}} \in \mathbb {R}^{{N_d}\times 4}$$ is built. Its row number refers to the global index of a data point (i). In each row, there are four elements, in turn representing the global index of the frequency ($${\mathcal {I}_{i1}} \in [1,N_f]$$), transmitter ($${\mathcal {I}_{i2}} \in [1,N_s]$$), receiver ($${\mathcal {I}_{i3}} \in [1, N_r]$$) and data type ($${\mathcal {I}_{i4}}$$). To calculate the ith data point $$d^\mathrm{syn}_i(\mathbf {m}^{})$$ on the current parameter model m, the real physical model $${\boldsymbol{\sigma }^{}}$$ must be obtained via the inverse model transformation (8b) first. Second, the discrete forward problem (5) with the model $${\boldsymbol{\sigma }^{}}$$ due to the frequency-transmitter pair $$(\omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})$$ is solved for the solution $${{\mathbf {x}}}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})$$. Finally, $$d^\mathrm{syn}_i(\mathbf {m}^{})$$ can be obtained by combing the primary field computed analytically in the background model with the interpolated finite-element solution   \begin{eqnarray} d^\mathrm{syn}_i(\mathbf {m}^{}) &=& d^p_i(\omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}}, {\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, \Theta _{{\mathcal {I}_{i4}}}) \nonumber\\ &&+ \mathbf {q}^T(\omega _{\mathcal {I}_{i1}}\!,{\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, \Theta _{{\mathcal {I}_{i4}}}) {{\mathbf {x}}}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}}), \end{eqnarray} (11)where $$\mathbf {q}(\omega _{\mathcal {I}_{i1}}\!,{\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, \Theta_{{\mathcal {I}_{i4}}}) \in \mathbb {C}^{N_{\mathrm{dof}}}$$ is the (complex) interpolation vector built to evaluate the field component $$\Theta_{{\mathcal {I}_{i4}}}$$ at the frequency $$\omega _{\mathcal {I}_{i1}}$$ and receiver $${\mathbf {r}^r_{{\mathcal {I}_{i3}}}}$$. It does not depend on the transmitter position and is constructed according to the interpolation rules of the curl-conforming finite elements and thus includes only geometric weights and frequency-related values for specific data types. For instance, if the datum is the x-component of the electric field, the specific expression of q reads   $$\mathbf {q}(\omega _{\mathcal {I}_{i1}}\!,{\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, E_x) = \left\lbrace N_{1,x}({\mathbf {r}^r_{{\mathcal {I}_{i3}}}}),\ldots , N_{\mathrm{dof},x}({\mathbf {r}^r_{{\mathcal {I}_{i3}}}})\right\rbrace ^T.$$ (12)Because the basis functions have a small support, q turns out to be very sparse. 3.1 An inexact Gauss–Newton framework To find a minimizer of the objective functional (9), we employ an inexact Gauss–Newton approach. In terms of the data Jacobian, the gradient ∇Φ(m) and the Gauss–Newton Hessian HGN(m) of the functional can be expressed as   \begin{eqnarray} {\nabla \Phi (\mathbf {m}) = }\nonumber\\ &&\, && \Re \left\lbrace \begin{array}{c}-[\mathbf {W}{\mathbf {J}_h(\mathbf {m}^{})}]^H \mathbf {W}[\mathbf {d}^\mathrm{obs}- \mathbf {d}^{\mathrm{syn}}(\mathbf {m}^{})] + \lambda _h \mathbf {L}_h^T\mathbf {L}_h(\mathbf {m}^{}_h - \mathbf {m}^{\mathrm{ref}}_h) \\ \quad -[\mathbf {W}{\mathbf {J}_v(\mathbf {m}^{})}]^H \mathbf {W}[\mathbf {d}^\mathrm{obs}- \mathbf {d}^{\mathrm{syn}}(\mathbf {m}^{})] + \lambda _v \mathbf {L}_v^T\mathbf {L}_v(\mathbf {m}^{}_v - \mathbf {m}^{\mathrm{ref}}_v) \end{array} \right\rbrace , \nonumber\\ \end{eqnarray} (13a)  \begin{eqnarray} {{\mathbf {H}^{\rm GN}(\mathbf {m})}}\nonumber\\ && {{ \!=\! \Re \!\left\lbrace \begin{array}{c}\left[\mathbf {W}{\mathbf {J}_h(\mathbf {m}^{})}\right]^H [\mathbf {W}{\mathbf {J}_h(\mathbf {m}^{})}] \!+\! \lambda _h \mathbf {L}_h^T\mathbf {L}_h \quad\ [\mathbf {W}{\mathbf {J}_h(\mathbf {m}^{})}]^H [\mathbf {W}{\mathbf {J}_v(\mathbf {m}^{})}] \\ \quad \left[\mathbf {W}{\mathbf {J}_v(\mathbf {m}^{})}\right]^H [\mathbf {W}{\mathbf {J}_h(\mathbf {m}^{})}] \quad\ [\mathbf {W}{\mathbf {J}_v(\mathbf {m}^{})}]^H [\mathbf {W}{\mathbf {J}_v(\mathbf {m}^{})}]\!+\!\lambda _v \mathbf {L}_v^T\mathbf {L}_v \end{array} \right\rbrace ,}}\nonumber\\ \end{eqnarray} (13b) where $${\mathbf {J}_h(\mathbf {m}^{})} \in \mathbb {C}^{N_d\times N_m}$$ and $${\mathbf {J}_v(\mathbf {m}^{})} \in \mathbb {C}^{N_d \times N_m}$$ are the data Jacobian components with respect to $$\mathbf {m}^{}_h$$ and $$\mathbf {m}^{}_v$$, respectively,   \begin{eqnarray} \left. {\mathbf {J}_h(\mathbf {m}^{})} \right|_{ij} &=& \frac{\partial d^\mathrm{syn}_i(\mathbf {m}^{})}{\partial {m}^{}_{h_{j}}}=-\mathbf {q}^T(\omega _{\mathcal {I}_{i1}}\!,{\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, \Theta _{{\mathcal {I}_{i4}}}) \mathbf {K}^{-1}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}) \nonumber \\ \quad \times\left[\frac{\partial \mathbf {K}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}})}{\partial {\sigma }^{}_{h_{j}}} {{\mathbf {x}}}({\boldsymbol{\sigma }},\omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})\right.\nonumber\\ \quad \left. - \frac{\partial {\mathbf {F}^{{}}}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})}{\partial {\sigma }^{}_{h_{j}}}\right]{\frac{\partial {\sigma }^{}_{h_{j}}}{\partial {m}^{}_{h_{j}}}}, \end{eqnarray} (14a)  \begin{eqnarray} \left. {\mathbf {J}_v(\mathbf {m}^{})} \right|_{ij} &=& \frac{\partial d^\mathrm{syn}_i(\mathbf {m}^{})}{\partial {m}^{}_{v_{j}}} = -\mathbf {q}^T(\omega _{\mathcal {I}_{i1}}\!,{\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, \Theta _{{\mathcal {I}_{i4}}}) \mathbf {K}^{-1}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}) \nonumber \\ \quad \times \left[\frac{\partial \mathbf {K}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}})}{\partial {\sigma }^{}_{v_{j}}} {{\mathbf {x}}}({\boldsymbol{\sigma }},\omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})\right.\nonumber\\ \quad -\left. \frac{\partial {\mathbf {F}^{{}}}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})}{\partial {\sigma }^{}_{v_{j}}}\right]{\frac{\partial {\sigma }^{}_{v_{j}}}{\partial {m}^{}_{v_{j}}}}. \end{eqnarray} (14b) We note that the mixed partial derivatives of the stabilizing functionals vanish because the regularizations on $$\mathbf {m}^{}_h$$ and $$\mathbf {m}^{}_v$$ are decoupled. Moreover, we see that the off-diagonal matrices of the Hessian (13b) actually couple the solutions through anisotropy. The Gauss–Newton method minimizes the non-linear functional (9) through solving the following linearized subproblem iteratively   $$\mathbf {H}^{\rm GN}(\mathbf {m}^{k}) \left( \begin{array}{c}\delta \mathbf {m}^{k}_h\\ \delta \mathbf {m}^{k}_v \end{array} \right) = -\nabla \Phi (\mathbf {m}^{k}),$$ (15)where k denotes the kth Gauss–Newton iteration. For practical 3-D CSEM inverse problems, this dense equation system has to be solved iteratively which only requires matrix-vector operations. To distinguish the approach for solving eq. (15) from the one for eq. (5), we denote the former by the inner solver. The normal eq. (15) might become very ill-conditioned for large-scale inverse problems with a non-uniform parametrization. Rather than resorting to the conventional CG method, a numerically more stable strategy is to solve the system (15) in the linear least-squares sense using some least-squares solver such as LSQR (Paige & Saunders 1982) as it only operates on J instead of JHJ  \begin{eqnarray} {\left( \begin{array}{c@\quad c}\Re \lbrace \mathbf {W}{\mathbf {J}_h(\mathbf {m}^{k})}\rbrace & \Re \lbrace \mathbf {W}{\mathbf {J}_v(\mathbf {m}^{k})}\rbrace \\ \Im \lbrace \mathbf {W}{\mathbf {J}_h(\mathbf {m}^{k})}\rbrace & \Im \lbrace \mathbf {W}{\mathbf {J}_v(\mathbf {m}^{k})}\rbrace \\ \sqrt{\lambda _h^k} \mathbf {L}_{h} & \mathbf {0} \\ \mathbf {0} & \sqrt{\lambda _v^k} \mathbf {L}_{v}\\ \end{array} \right) \left( \begin{array}{c}\delta \mathbf {m}^{k}_h\\ \delta \mathbf {m}^{k}_v \end{array} \right)} \nonumber\\ =\left( \begin{array}{cc}\Re \lbrace \mathbf {W}[\mathbf {d}^{\mathrm{obs}}- \mathbf {d}^{\mathrm{syn}}(\mathbf {m}^{k})]\rbrace \\ \Im \lbrace \mathbf {W}[\mathbf {d}^{\mathrm{obs}}- \mathbf {d}^{\mathrm{syn}}(\mathbf {m}^{k})]\rbrace \\ -\sqrt{\lambda _h^k} \mathbf {L}_{h}(\mathbf {m}^{k}_h-\mathbf {m}^{\mathrm{ref}}_h) \\ -\sqrt{\lambda _v^k} \mathbf {L}_{v}(\mathbf {m}^{k}_v-\mathbf {m}^{\mathrm{ref}}_v) \end{array}\right). \end{eqnarray} (16) For large-scale problems, on the one hand, solving eq. (16) exactly is impractical. On the other hand, it is rather unnecessary to obtain an exact solution especially when the current model mk is far away from a (local) minimum. Some studies suggest that regularization is implicitly applied when the solving procedure is terminated at early iterations (Santos 1996; Kilmer & O’Leary 2001; Grayver et al.2013). The truncation of the iteration is problem-dependent and many authors have found some useful stopping criteria by empirical investigations to explain the observations to the level of measurement errors (e.g. Grayver et al.2013; Schwarzbach & Haber 2013). The Gauss–Newton method with an early termination when solving the linearized subproblem (16) is termed the inexact Gauss–Newton method (Nocedal & Wright 1999). To prevent the Gauss–Newton method from being divergent, an inexact line search algorithm according to the Armijo condition is augmented to steer the decrease in the objective functional (Nocedal & Wright 1999). 3.2 The data Jacobian In a Gauss–Newton-based inversion scheme, preparing the full data Jacobian at each iteration is the most time-consuming ingredient besides solving the forward problem. Determining an efficient strategy for treating the Jacobian is problem-dependent. To accommodate most of the cases, we have provided three alternatives in this framework to deal with the Jacobian. For small- or moderate-size problems, the Jacobian can be constructed explicitly. According to eq. (14), this will lead to a memory cost of $$\mathcal {O}(N_d\times N_m\times 2)$$ by carrying out Nf × Nr × Nt additional forward runs. With an explicit Jacobian, operating the Jacobian-vector products becomes very efficient which allows the inner solver to carry out more iterations for a more accurate solution (Scheunert et al.2016). In contrast to the explicit storage, a widely used strategy is to treat the Jacobian implicitly which works for an iterative solution of the Gauss–Newton equation system (15). If the inner iterative process is terminated after Nitr iterations, 2 × Nf × Ns × Nitr additional forward runs are needed to obtain a model update δmk. This approach is memory-efficient, but time-consuming because the required number of forward simulations increases linearly with the number of inner iterations Nitr. Besides the above two approaches, we introduce a third method which we call the semi-explicit strategy. It acts as a compromise between the explicit and implicit methods because only one part of the Jacobian is stored explicitly which requires forward runs. At each Gauss–Newton iteration, we first compute and store the vectors qTK − 1 occurring in the Jacobian (i.e. the solution of the adjoint problem at each receiver, see eq. 14) by exploiting the symmetry of the system matrix K, which needs a memory of $$\mathcal {O}(N_f \times N_r \times N_t \times N_{\mathrm{dof}})$$ and, like the explicit storage, Nf × Nr × Nt additional forward runs (Newman & Alumbaugh 1997). Any Jacobian-vector multiplication Jp can be done by first carrying out a sparse matrix-vector product $$(\frac{\partial \mathbf {K}}{\partial \boldsymbol{\sigma }} \mathbf {x} - \frac{\partial \mathbf {F}}{\partial \boldsymbol{\sigma }})\frac{\partial \boldsymbol{\sigma }}{\partial \mathbf {m}} \mathbf {p}$$ and then multiplying the pre-computed row vectors by the resulting vector. The inverse procedure applies when operating the transpose of the data Jacobian. The memory cost in the semi-explicit method might be less demanding than that in the explicit storage if there are only tens to hundreds of receiver positions which can be the case in marine CSEM applications (Gabrielsen et al.2009). The efficiency of Jacobian-vector operations compared to that in the explicit method largely relies on the relation between Nm × 2 and Ndof. For high-order curl-conforming basis functions, the latter can be much larger than the former one, while they might be comparable when the lowest order basis function is considered. Nonetheless, a notable advantage of this strategy is that the memory requirement is not affected by the number of model parameters. This makes it competitive for large-scale anisotropic 3-D CSEM inverse problems in which the number of model parameters increases linearly with the dimension of anisotropy. 4 TEST AND APPLICATION In this section, we first examine the new anisotropic inversion scheme using a synthetic marine model with rough bathymetry. Based on this, we apply it to interpret the Troll field data acquired in 2008 by Statoil ASA and EMGS ASA in the North Sea (Gabrielsen et al.2009). All computations were conducted on a moderately parallel machine equipped with 4 Intel(R) Xeon(R) 2.2 GHz CPUs with 8 cores per CPU and 512 GB of RAM. To accelerate the inversion process, we have implemented a parallel scheme typically designed for marine CSEM data at a few frequencies and many transmitter positions. First, the computational model is built and discretized considering all frequencies, transmitters and receivers. Then, the entire data set is decomposed with respect to frequencies. As solving the system of linear equations for the forward problem with a large number of transmitters can be a significant burden when using a direct solver, each data subset can be further decomposed with respect to transmitters to accelerate the solution phase. Each computational task which consists of computing the Jacobian-vector products, the gradient and misfit of the objective function for the corresponding data subset is then run in parallel using the Message Passing Interface specification. Inside each process, the finite-element equation systems are solved using MUMPS (Amestoy et al.2001). For a direct solver, the decomposition with respect to transmitters is optional when sufficient memory is available. For the synthetic example, we decompose the data set with respect to frequencies. For the field data set, the parallel decomposition is also performed with respect to the transmitters. The inner solver LSQR is terminated when either the maximum number of iterations (120 for the synthetic example and 60 for the field experiment) is achieved or the tolerance 10−2 is satisfied. The regularization parameters λh and λv are chosen to be equal unless otherwise stated. We terminate the inversion process once the data misfit does not decrease significantly. 4.1 A synthetic marine CSEM model The model consists of an air layer of 10−7 Sm−1, a sea water column of 3.7037 Sm−1, an undulating seafloor, a subsurface with a horizontal conductivity of 1 Sm−1 and a vertical conductivity of 0.25 Sm−1, and two resistive anomalies of 0.01 Sm−1. Considering a right-handed Cartesian coordinate system with the z-axis pointing downwards, the flat sea surface is situated at z = 0 m. The depth of the sea water in the survey area varies from about z = 352 m to about z = 916 m (see the depth contour in Fig. 1a). The smaller resistive body has a dimension of 1500 m × 2000 m × 100 m along the x-, y- and z-directions with the centre being located at (2750 m, −500 m, 1400 m). The larger body has a size of 3000 m × 3000 m × 100 m and is centred at (− 1500 m, 0 m, 1500 m). Figure 1. View largeDownload slide Survey geometry over the undulating seafloor. (a) It displays the depth contour of the seafloor topography and the dashed black lines outline the horizontal extents of the anomalies. Receivers are denoted by black triangles and transmitters by red dots. (b) It depicts the vertical cross-section of the subsurface along transmitter line Tx03 showing the depth extents of the two anomalous bodies. (c) It highlights the bathymetry variation along each towline Tx01 to Tx04. Figure 1. View largeDownload slide Survey geometry over the undulating seafloor. (a) It displays the depth contour of the seafloor topography and the dashed black lines outline the horizontal extents of the anomalies. Receivers are denoted by black triangles and transmitters by red dots. (b) It depicts the vertical cross-section of the subsurface along transmitter line Tx03 showing the depth extents of the two anomalous bodies. (c) It highlights the bathymetry variation along each towline Tx01 to Tx04. We consider a dense CSEM survey including four transmitter lines named Tx01 to Tx04. As shown in Fig. 1, 44 seabed receivers (black triangles) are arranged in an 11 × 4 rectangular array with an x-directed spacing of 900 m and a y-directed spacing of 1000 m. Each receiver is equipped with two orthogonal electric antennae for recording the horizontal components of the electric field $$E_{x^{\prime }}$$ and $$E_{y^{\prime }}$$. Due to the bathymetry, the sensors generally have a non-zero pitch (right-handed rotation around the y-axis, positive bow up) which varies from −15.8° to 12.0° as shown in Fig. 2. Consequently, $$E_{x^{\prime }}$$ and $$E_{y^{\prime }}$$ might not be horizontally oriented. We do not rotate the components into the global coordinate system but denote and compute them in each receiver’s local coordinates using the exact rotation angles. The vertical position of the receivers varies from around z = 438 m to around z = 794 m. An x-directed HED of 200 m length is towed about 30 m above each x-directed receiver line ranging from x = −8100 to 8100 m with an acquisition separation of 300 m (red dots in Fig. 1). This gives rise to a total number of 220 transmitter positions. The transmission frequencies were 0.25 and 0.75 Hz. Both inline and offline data within the transmitter–receiver range between 500 and 9000 m were recorded yielding 33 792 complex measurements in total. Figure 2. View largeDownload slide Receiver pitch. Figure 2. View largeDownload slide Receiver pitch. Since the number of transmitters (220) is larger than the product of the number of receivers (44) and the number of data types (2), we apply the reciprocity principle (Zhdanov 2009). By doing this, the real receiver antennae will become the computational transmitters. Correspondingly, the computational receivers are now placed where the real transmitters are located to record the component along the transmitter’s orientation. This leads to 44 x΄-directed and 44 y΄-directed computational electric transmitters (both in the local coordinate system) and 220 x-directed electric receivers with a length of 200 m. We first prepare the model for generating the synthetic data set. It contains the two anomalies embedded in the rugged subsurface with a size of about 60000 m × 60000 m × 15000 m, the sea water column and an air layer with a thickness of 70000 m on top. The whole computational domain is decomposed into 746662 tetrahedra. The second-order finite-element method is employed to solve the forward problem yielding 4724736 degrees of freedom. The background model is chosen to be the air–water half-space (the same has been used for the inversion). Although not shown here, we have checked the modeling accuracy on this mesh by comparing to the simulation result of a planar layer. The synthetic data were normalized by the dipole moment (200 Am) and 3 per cent Gaussian noise was added. To avoid the inverse crime, we set up a new mesh for the inversion. The size of the computational domain stays the same but the geometries of the two anomalies are excluded (Fig. 3). Since the success of a pixel-based parametrization strategy greatly relies on a proper discretization of the inversion domain, we refine the subsurface below the receiver grid uniformly from about z = 1100 to 1700 m to ensure a sufficient resolution of the thin resistors. The computational domain is discretized into 672 154 tetrahedra. The forward operator is based on the second-order finite-element discretization giving 4255 588 degrees of freedom. We consider the entire subsurface which includes 427 509 cells as the inversion domain. We set the lower and upper conductivity bounds to 0.002 and 10 Sm−1, respectively, in the model transformation (8). For this setup, the explicit storage of the anisotropic Jacobian requires 33 792 × 427509 × 2 × 16 B ≈ 430.5 GB of memory which is not affordable on the current computing system. In contrast, the semi-explicit storage strategy merely needs 2 × 220 × 4255588 × 16 B ≈ 27.9 GB of memory. Table 1 lists a comparison of all three methods for treating the Jacobian in terms of the run times and memory cost. From this example, it can be seen that the semi-explicit storage outperforms the implicit storage in terms of speed and the explicit storage in terms of memory requirement. Figure 3. View largeDownload slide Model discretization for synthetic data generation (upper panel) and the inversion (lower panel) at (a) and (d) x = 0 m, (b) and (e) y = 0 m and (c) and (f) z = 1400 m. Figure 3. View largeDownload slide Model discretization for synthetic data generation (upper panel) and the inversion (lower panel) at (a) and (d) x = 0 m, (b) and (e) y = 0 m and (c) and (f) z = 1400 m. Table 1. The performance of different strategies for treating the Jacobian.   Explicit  Implicit  Semi-explicit  One Jacobian-vector product (s)  –  170.5  13.5  Storage (GB)  430.5   0   27.9     Explicit  Implicit  Semi-explicit  One Jacobian-vector product (s)  –  170.5  13.5  Storage (GB)  430.5   0   27.9   View Large Moreover, it is worthwhile to mention that, in the case of the semi-explicit storage, the memory cost for storing the forward solution vectors and the vectors qTK − 1 is the same as used without applying the reciprocity principle. However, the use accelerates the computation of the primary fields and the assembly of the source terms of the forward problem (88 electric dipoles with reciprocity in comparison with 220 electric line sources without reciprocity). The first inversion experiment (E0) starts with a homogeneous isotropic subsurface of 0.5 Sm−1 which differs from the true anisotropic background (σh = 1 Sm−1 and σv = 0.25 Sm−1). We choose the reference model to be the same as the starting model (also for the following four tests). The initial value of the regularization parameters $$\lambda _h^0$$ and $$\lambda _v^0$$ is 0.000164. We decrease them by a factor of 0.6 at each iteration. Fig. 4 shows the slices of the final inversion results. For the sake of illustration, all inversion results are displayed in terms of the electrical resistivity. We observe from this example that a comprehensive data coverage with both inline and azimuthal data can determine the electrical anisotropy in the sediment. The maximum value of the resolved vertical resistivity in the reservoirs is about 30 Ωm compared to the true value 100 Ωm. But the reconstructed resistors are also thicker as seen in the vertical sections. Therefore, it results in a reasonable anomalous transverse resistance (about 11495 Ωm2 for the large body and 4389 Ωm2 for the small one) compared to the true value (9600 Ωm2 for both targets). Anomalous transverse resistance is the quantity that determines the CSEM response. It is also seen from Fig. 4 that the model was updated in the regions with weak sensitivity which is the effect of the smoothness operator. Figs 5 and 6 display the relative amplitude and phase residuals in terms of the data uncertainty. The comparison with the ‘exact’ relative data residuals for the noise-free synthetic data [denoted by (d*, ϕ*)] indicates that the predicted model fits the observed data well in the order of magnitude of the noise. Figure 4. View largeDownload slide Slices of the resolved horizontal (left-hand panel) and vertical (right-hand panel) resistivity models in the experiment E0 at (a) and (b) z = 1480 m, (c) and (d) z = 1360 m, along towline (e) and (f) Tx01, (g) and (h) Tx02, (i) and (j) Tx03 and (k) and (l) Tx04. Figure 4. View largeDownload slide Slices of the resolved horizontal (left-hand panel) and vertical (right-hand panel) resistivity models in the experiment E0 at (a) and (b) z = 1480 m, (c) and (d) z = 1360 m, along towline (e) and (f) Tx01, (g) and (h) Tx02, (i) and (j) Tx03 and (k) and (l) Tx04. Figure 5. View largeDownload slide Relative amplitude misfit in the experiment E0 at (a) the 0th iteration and (b) the 13th iteration in comparison with (c) the true misfit for noise-free data (d*, ϕ*). Transmitters and receivers are indexed by firstly walking along the positive y-direction and then the positive x-direction. Figure 5. View largeDownload slide Relative amplitude misfit in the experiment E0 at (a) the 0th iteration and (b) the 13th iteration in comparison with (c) the true misfit for noise-free data (d*, ϕ*). Transmitters and receivers are indexed by firstly walking along the positive y-direction and then the positive x-direction. Figure 6. View largeDownload slide Relative phase misfit in the experiment E0 at (a)  the 0th iteration and (b)  the 13th iteration in comparison with (c) the true misfit for noise-free data (d*, ϕ*). Figure 6. View largeDownload slide Relative phase misfit in the experiment E0 at (a)  the 0th iteration and (b)  the 13th iteration in comparison with (c) the true misfit for noise-free data (d*, ϕ*). Figure 7. View largeDownload slide The development of the (a) data misfit function and (b)  stabilizing functions for different inversion conditions in the synthetic model studies. Figure 7. View largeDownload slide The development of the (a) data misfit function and (b)  stabilizing functions for different inversion conditions in the synthetic model studies. To further investigate the dependence of the inversion results on the starting model, the number of inner iterations Nitr, and the model constraints, we carry out another four experiments. To examine whether the current inversion result is affected by insufficiently solving each Gauss–Newton system of equations (16) or a loose model constraint, we perform the first additional model study (E1) by setting Nitr = 150 and a stricter model bound of σ ∈ (0.005, 2) Sm−1. The second (E2) and third (E3) additional experiments examine the influence of different starting models. E2 considers a homogeneous subsurface of σh = σv = 1 Sm−1 as the starting model with Nitr = 120 and a bound constraint σ ∈ (0.002, 10) Sm−1. E3 utilizes the same inversion conditions for E2 except for the true background model as the initial guess. The fourth additional experiment (E4) is based on experiment E3 but using Nitr = 300. Figs 7a and 7b plot the development of the data misfit function Φd and regularization functions $$\Phi _{m_h}$$ and $$\Phi _{m_v}$$ for the above five experiments, respectively. Despite the use of different inversion conditions, the five tests have converged to the same level: E0 converges to a total data misfit of 0.5520, E1 is 0.5531, E2 is 0.5493, E3 is 0.5515 and E4 is 0.5511. Fig. 8 depicts the resistivities along vertical pseudo-boreholes as a function of depth at the horizontal receiver locations of (−1800 m, 0 m) and (2700 m, −1000 m). We see that a large number of LSQR iterations in E4 does not significantly improve the inversion result indicating Nitr = 120 is sufficient for the current example. The stronger model constraint used in E1 seems to enhance the resolved anomalies compared to other setups. Figure 8. View largeDownload slide Resistivities along a vertical borehole through the anisotropic inversion models at the horizontal positions of (a) and (b) (2700 m, −1000 m) and (c) and (d)( − 1800 m, 0 m). Figure 8. View largeDownload slide Resistivities along a vertical borehole through the anisotropic inversion models at the horizontal positions of (a) and (b) (2700 m, −1000 m) and (c) and (d)( − 1800 m, 0 m). We also note that the horizontal resistivity of the two anomalies (100 Ωm) is hardly resolved with this survey configuration. This phenomenon is caused by the guided wave excited in the thin resistive layer. Due to the continuity equation, a very strong vertical electric field is generated inside the resistor. The magnitude of the horizontal electric field is in comparison negligible which prevents the horizontal resistivity from being resolved (Ramananjaona et al.2011; Mittet & Morten 2013). Moreover, Brown et al. (2012) pointed out that using an HED source, the horizontal resistivity of a thin resistive reservoir cannot be reconstructed once its value is greater than that of the background. This finding corresponds well with the observation here. 4.2 Troll field The Troll field is located in a mature area in the North Sea where early CSEM surveys showed that the CSEM method can be a powerful exploration technique due to the prominent variation in EM fields caused by the large-scale resistive gas reservoir. Some authors have reported the imaging results of the CSEM data set acquired in 2003 using a single receiver line across the Troll West Gas Province (TWGP, e.g. Commer & Newman 2008; Plessix & Mulder 2008). Here, we utilize a more comprehensive CSEM data set collected in 2008 using a dense 3-D acquisition grid by Statoil ASA and EMGS ASA to map the resistivity distribution of this area (Gabrielsen et al.2009; Morten et al.2012; Jaysaval et al.2014). The whole survey contains six north–south towlines being around 1.25 km apart, three east–west towlines with a large line spacing of about 3.8 km, and 53 stationary seabed detectors aiming to focus on the relatively small-scale Troll West Oil Province (TWOP) while covering only a small part of the TWGP (see Gabrielsen et al.2009). Both inline and azimuthal data were recorded simultaneously and processed at frequencies 0.25, 0.75 and 1.25 Hz. The inline measurements from the westernmost receiver line were excluded because of the large impact of pipelines. For a clear illustration, we only show the processed layout in Fig. 9 by excluding the westernmost receiver line. The processed data at hand are the rotated horizontal components of the electric field Ex and Ey with Ex being along the corresponding towline direction and Ey being perpendicular to that. The data uncertainties of the rotated field components are related to the data uncertainties of the non-rotated field components through error propagation analysis and the data weights are constructed based on that according to Morten et al. (2009). Figure 9. View largeDownload slide The adapted Troll field survey layout in 2008 (Gabrielsen et al.2009). The colour map describes the variation of the water depth. Black rectangles indicate seabed receivers and red lines represent towlines. The TWOP is outlined in white and TWGP in black. Figure 9. View largeDownload slide The adapted Troll field survey layout in 2008 (Gabrielsen et al.2009). The colour map describes the variation of the water depth. Black rectangles indicate seabed receivers and red lines represent towlines. The TWOP is outlined in white and TWGP in black. For the inversion, we select the inline and azimuthal data from the five north–south and three east–west towlines at frequencies 0.25 and 0.75 Hz within the transmitter–receiver range between 1200 and 9000 m. The data subset includes 157 748 data from 44 Ex and 44 Ey receivers and 38 349 transmitter positions. Again, we apply the reciprocity principle resulting in 88 computational seabed transmitters and 38 349 towline-directed electric receivers of 270.5 m length. The bathymetry over the survey area (Fig. 9) is rather smooth and small compared to the survey scale. The water depth varies from about 300 m in the southeast to about 360 m in the northeast. Although previous studies have assumed that the bathymetry over this region does not influence the measurements much (Newman et al.2010), we incorporated it into the model to highlight the flexibility of unstructured meshes and the finite-element approach. We build and visualize the inversion models using a coordinate system with x along the north direction, y along the east direction and z pointing downwards. The sea surface is at z = 0 m. The subsurface with a topographic seafloor has a size of about 70000 m × 66000 m × 30000 m, overlain by sea water with a conductivity of 3.7037 Sm−1, and an air layer of 10−7 Sm−1. We prepared distinct meshes for isotropic and anisotropic inversions, respectively. This is reasonable because our preliminary study showed that many small-scale structures were generated close to the seafloor in the isotropic inversion. Therefore, it shall require a denser discretization in that case to rule out the possibility that this is caused by the inaccuracy of the forward simulation. In order to ensure a reasonable resolution for both oil and gas fields, we refine the subdomain which laterally covers the horizontal extents of the oil and gas fields and vertically extends from the sea bottom to about z = 2700 m. The vertical spacing of the cells in the refined region is around 100 m. Outside the subdomain, the mesh gets coarser towards the outer boundaries. The regions around the transmitters and receivers are locally refined with a minimum edge length of around 0.4 m to ensure sufficient accuracy. For anisotropic inversions, the whole computational domain is decomposed into 4882 262 tetrahedra, 2269 297 of which are distributed in the sea water layer and 2337 595 in the subsurface. The cells close to the seafloor are further refined in the isotropic case. The resulting denser mesh consists of 5152 631 cells, 2401 811 of which are in the sea water and 2463 768 in the subsurface. For all experiments, we employ the first-order finite-element method to solve the forward problem with an air–water half-space as the background model. The whole subsurface is considered as the inversion domain. The model constraint is set to σ ∈ (0.002, 5) Sm−1. Other inversion setups and statistics are listed in Table 2. For this large-scale configuration, neither the explicit (too many model parameters) nor the semi-explicit (too many computational receivers) storage of the data Jacobian is feasible on our current computing system. Therefore, the Jacobian is treated implicitly. Each inversion task is decomposed into 2 × 2 (number of partitions over frequencies × number of partitions over transmitters) inversion jobs. With a direct solver, the drawback of this decomposition is the repeated storage of the factor of the system matrix at the same frequency. However, if the memory requirement is permitted, a scalability factor of (nearly) 2 is observed compared to a 2 × 1 decomposition. Table 2. Inversion conditions for the three Troll experiments.   Isotropic  Anisotropic (inline)  Anisotropic (inline and azimuthal)  Mesh size  5152 631  4882 262  4882 262  Nm  2463 768  2337 595  2337 595  Nd  16 088  16 088  157 748  Ndof  5985 264  5672 052  5672 052  $$\lambda _h^0$$ and $$\lambda _v^0$$  0.000252  0.000207  0.0000255  mrefh and mrefv  σ = 0.5 Sm−1  σh = σv = 0.5 Sm−1  σh = σv = 0.5 Sm−1  Run time (hour)  59.1  48.2  67.4  Matrix factor (GB)  67.3  62.6  62.6    Isotropic  Anisotropic (inline)  Anisotropic (inline and azimuthal)  Mesh size  5152 631  4882 262  4882 262  Nm  2463 768  2337 595  2337 595  Nd  16 088  16 088  157 748  Ndof  5985 264  5672 052  5672 052  $$\lambda _h^0$$ and $$\lambda _v^0$$  0.000252  0.000207  0.0000255  mrefh and mrefv  σ = 0.5 Sm−1  σh = σv = 0.5 Sm−1  σh = σv = 0.5 Sm−1  Run time (hour)  59.1  48.2  67.4  Matrix factor (GB)  67.3  62.6  62.6  View Large First, we carry out an isotropic inversion of the inline data from the selected data subset. The left-hand panel of Fig. 10 shows the slices of the resolved resistivity model at the depth of z = 1500 m and along the five north–south towlines. Although both oil and gas fields are indicated, there exist many small-scale anomalies which likely are imaging artefacts if compared with a priori geostructures. This poses the question of whether electrical anisotropy is present. To demonstrate this assumption, secondly, we run an anisotropic inversion of the same inline data subset and display the recovered vertical resistivity model in the right-hand column of Fig. 10. The comparison shows that the anisotropic inversion renders much fewer artefacts around the receivers and near the seafloor than the isotropic inversion. Moreover, the anisotropic resistivity model becomes geologically more plausible and the horizontal extents of the resistivity anomalies coincide well with the outlined a priori regions. The oscillating layers observed in the overburden in the isotropic inversion, which can to some extent represent an effect like anisotropy, have largely disappeared from the anisotropic inversion. To show an example how the models fit the measurements, the relative data residuals at the common offset of 2500 m are checked for both inversions. They are positioned at the common midpoint of the corresponding transmitters and receivers. It is seen from Figs 11 and 12 that the anisotropic resistivity model yields a better data fit than the isotropic resistivity model. Based on the above observations, electrical anisotropy seems to be significant in this data set. Figure 10. View largeDownload slide Slices of the isotropic (left-hand panel) and anisotropic (right-hand panel) Troll resistivity models characterized by merely inline data at (a) and (b)  z = 1500 m, along towline (c) and (d)  Tx002, (e) and (f) Tx003, (g) and (h) Tx004, (i) and (j) Tx005 and (k) and (l) Tx006. The solid black and white lines in (a) and (b) outline the horizontal extents of the Troll gas and oil fields, respectively, and black triangles indicate receiver positions. Figure 10. View largeDownload slide Slices of the isotropic (left-hand panel) and anisotropic (right-hand panel) Troll resistivity models characterized by merely inline data at (a) and (b)  z = 1500 m, along towline (c) and (d)  Tx002, (e) and (f) Tx003, (g) and (h) Tx004, (i) and (j) Tx005 and (k) and (l) Tx006. The solid black and white lines in (a) and (b) outline the horizontal extents of the Troll gas and oil fields, respectively, and black triangles indicate receiver positions. Figure 11. View largeDownload slide Relative amplitude misfit from the isotropic (left-hand panel) and anisotropic (right-hand panel) inversions of the Troll inline data in the common-midpoint domain at the offset of 2500 m. Figure 11. View largeDownload slide Relative amplitude misfit from the isotropic (left-hand panel) and anisotropic (right-hand panel) inversions of the Troll inline data in the common-midpoint domain at the offset of 2500 m. Figure 12. View largeDownload slide Relative phase misfit from the isotropic (left-hand panel) and anisotropic (right-hand panel) inversions of the Troll inline data in the common-midpoint domain at the offset of 2500 m. Figure 12. View largeDownload slide Relative phase misfit from the isotropic (left-hand panel) and anisotropic (right-hand panel) inversions of the Troll inline data in the common-midpoint domain at the offset of 2500 m. Finally, we run another anisotropic inversion by additionally incorporating the azimuthal data. The slices of the imaged horizontal and vertical resistivity models are displayed in Fig. 13 and the selected data residuals in Fig. 14. It is seen through comparing Figs 13 and 10 that, despite the fact that the anisotropic inversion of purely inline data has yielded the main characteristics of the Troll oil and gas reservoirs, the incorporation of the azimuthal data has noticeably enhanced the reconstruction of both fields: (i) both the TWOP and TWGP are now indicated more clearly; (ii) their lateral extents are enlarged and (iii) the up-tilted part of the gas reservoir, which is strongly observed in the isotropic inversion while slightly shown in the corresponding anisotropic case, is now pushed down. According to the geological description and the inversion results shown by Gabrielsen et al. (2009) and Morten et al. (2012), the oil and gas fields have been resolved at the correct depth. The depth extents of both oil and gas fields are a bit too large ranging from about z = 1250 to about 1650 m for the oil field and from about z = 1100 to about 1800 m for the gas field. This is probably caused by the limited resolution of the data, the use of the smoothness regularization, and the relatively large vertical spacing of the cells (about 100 m). Furthermore, as shown in Fig. 15, the two anisotropic inversion processes steadily converge to a smaller overall data misfit than the isotropic inversion procedure. Figure 13. View largeDownload slide Slices of the anisotropic Troll resistivity model characterized by both inline and azimuthal data at depth (a) and (b) z = 1500 m, along towline (c) and (d) Tx002, (e) and (f) Tx003, (g) and (h) Tx004, (i) and (j) Tx005 and (k) and (l) Tx006. Figure 13. View largeDownload slide Slices of the anisotropic Troll resistivity model characterized by both inline and azimuthal data at depth (a) and (b) z = 1500 m, along towline (c) and (d) Tx002, (e) and (f) Tx003, (g) and (h) Tx004, (i) and (j) Tx005 and (k) and (l) Tx006. Figure 14. View largeDownload slide Relative amplitude (left-hand panel) and phase (right-hand panel) misfits from the anisotropic inversion of the full Troll data subset in the common-midpoint domain at the offset of 2500 m. Figure 14. View largeDownload slide Relative amplitude (left-hand panel) and phase (right-hand panel) misfits from the anisotropic inversion of the full Troll data subset in the common-midpoint domain at the offset of 2500 m. Figure 15. View largeDownload slide The development of the Troll (a) data and (b)  regularization misfits at each Gauss–Newton iteration. Note that the curve for $$\Phi _{m_h}$$ is absent in the isotropic case. Figure 15. View largeDownload slide The development of the Troll (a) data and (b)  regularization misfits at each Gauss–Newton iteration. Note that the curve for $$\Phi _{m_h}$$ is absent in the isotropic case. 5 CONCLUSIONS In this paper, we have presented an anisotropic 3-D inversion framework based on a finite-element discretization on unstructured grids for interpreting CSEM data in the frequency domain. The complexity of finite-element modeling compared to that of finite-difference modeling is balanced by its flexibility in incorporating different types of meshes. The freedom that comes with an unstructured mesh allows to describe very complicated geometries and avoids interpolation problems that are typical for rectilinear mesh approaches. This advantage becomes even more attractive when simulating large 3-D CSEM data sets with many transmitter positions involved in complex geometries. However, the application of finite-element methods is, to some extent, impeded by the effort of generating a quality mesh for a very complex model. To deal with this issue, some promising workflows have been reported in which complicated 3-D surfaces are represented in terms of non-uniform rational basis splines (Börner et al.2015; Zehner et al.2015). The resulting free surface can therefore be flexibly triangulated to avoid mesh degeneration caused by the description of a 3-D undulating surface through interpolating a fixed triangulation. The inversion scheme has incorporated two different formulations, namely by using the total and secondary field approaches for the forward problem, respectively, and offered three different strategies for treating the data Jacobian for a general application. Moreover, we have achieved a parallelism that allows us to tackle an industry-scale 3-D data set on a moderate-scale parallel machine. We restricted ourselves to marine CSEM applications in this paper. We picked the secondary field approach and validated the functionalities of the inversion scheme using a close-to-real-world synthetic marine CSEM scenario and the realistic Troll field data set from a 3-D dense survey grid in 2008. In this framework, we chose to use the parallel direct solver MUMPS for solving the system of linear equations for the multisource forward problem. Although direct solvers are oftentimes criticized because of high memory demand and poor scalability for large-scale problems, the Troll field example with about 5.7 million degrees of freedom has shown that the performance of the direct solver is still competitive. The memory requirement in a problem of such a scale is not only determined by the direct solver. The memory cost for storing the solutions, the sources (the right-hand side of the forward equation system), and some other auxiliary vectors for Jacobian-related operations is also non-trivial. However, when the problem size gets even larger, one might need to further optimize the forward solver by considering methods such as domain decomposition strategies (Rung-Arunwan & Siripunvaraporn 2010; Bakr et al.2013; Ren et al.2014) or/and fast convergent iterative approaches to relieve the memory request for storing the factor of the entire system matrix. We put this into our future research task list. Both the synthetic and real model studies have shown that the horizontal resistivity of thin resistive reservoirs cannot be resolved. The lack of sensitivity to the horizontal resistivity might be overcome by adding measurements due to other sources which generate a TE-mode dominant pattern in thin resistive anomalies. The significance of taking electrical anisotropy into account for marine CSEM data has been addressed by comparing the isotropic and anisotropic inversion results of the Troll field data. Acknowledgements We are grateful to EMGS ASA for providing the Troll field data set and for the permission to publish the results. FW would like to thank Christoph Schwarzbach for offering the finite-element forward operator. The work was partially funded by the China Scholarship Council (no. 2009637116). And finally, we thank editor Gary Egbert, the two anonymous reviewers, and Colin Farquharson for their insightful comments. REFERENCES Abubakar A., Liu J., Li M., Habashy T.M., MacLennan K., 2010. Sensitivity study of multi-sources receivers CSEM data for TI-anisotropy medium using 2.5D forward and inversion algorithm, in 72nd Annual International Conference and Exhibition , EAGE, Extended Abstracts, doi:10.3997/2214-4609.201400664 Amaya M., Morten J.P., Boman L., 2016. A low-rank approximation for large-scale 3D controlled-source electromagnetic Gauss-Newton inversion, Geophysics , 81( 3), E211– E225. Google Scholar CrossRef Search ADS   Amestoy P.R., Duff I.S., L’Excellent J.-Y., Koster J., 2001. A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM J. Matrix Anal. Appl. , 23( 1), 15– 41. Google Scholar CrossRef Search ADS   Anderson C., Mattson J., 2010. An integrated approach to marine electromagnetic surveying using a towed streamer and source, First Break , 28( 5), 71– 75. Ansari S., Farquharson C.G., 2014. 3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids, Geophysics , 79( 4), E149– E165. Google Scholar CrossRef Search ADS   Bakr S.A., Pardo D., Mannseth T., 2013. Domain decomposition Fourier finite element method for the simulation of 3D marine CSEM measurements, J. Comput. Phys. , 255, 456– 470. Google Scholar CrossRef Search ADS   Börner J.H., Wang F., Weißflog J., Bär M., Görz I., Spitzer K., 2015. Multi-method virtual electromagnetic experiments for developing suitable monitoring designs: A fictitious CO2 sequestration scenario in Northern Germany, Geophys. Prospect. , 63( 6), 1430– 1449. Google Scholar CrossRef Search ADS   Brown V., Hoversten M., Key K., Chen J., 2012. Resolution of reservoir scale electrical anisotropy from marine CSEM data, Geophysics , 77( 2), E147– E158. Google Scholar CrossRef Search ADS   Castillo P., Rieben R., White D., 2005. FEMSTER: an object-oriented class library of high-order discrete differential forms, ACM Trans. Math. Softw. , 31( 4), 425– 457. Google Scholar CrossRef Search ADS   Commer M., Newman G.A., 2008. New advances in three-dimensional controlled-source electromagnetic inversion, Geophys. J. Int. , 172( 2), 513– 535. Google Scholar CrossRef Search ADS   Commer M., Newman G.A., 2009. Three-dimensional controlled-source electromagnetic and magnetotelluric joint inversion, Geophys. J. Int. , 178( 3), 1305– 1316. Google Scholar CrossRef Search ADS   Constable S., 2010. Ten years of marine CSEM for hydrocarbon exploration, Geophysics , 75( 5), 75A67– 75A81. Google Scholar CrossRef Search ADS   Constable S., Srnka L.J., 2007. An introduction to marine controlled-source electromagnetic methods for hydrocarbon exploration, Geophysics , 72( 2), WA3– WA12. Google Scholar CrossRef Search ADS   Eidesmo T., Ellingsrud S., MacGregor L.M., Constable S., Sinha M.C., Johansen S.E., Kong F.N., Westerdahl H., 2002. Sea Bed Logging (SBL), a new method for remote and direct identification of hydrocarbon filled layers in deepwater areas, First Break , 20( 3), 144– 152. Ellis M.H., Sinha M., Parr R., 2010. Role of fine-scale layering and grain alignment in the electrical anisotropy of marine sediments, First Break , 28( 9), 49– 57. Fanavoll S., Gabrielsen P.T., Ellingsrud S., 2014a. The impact of CSEM on exploration decisions and seismic: two case studies from the Barents Sea, First Break , 32( 11), 105– 110. Fanavoll S., Gabrielsen P.T., Ellingsrud S., 2014b. CSEM as a tool for better exploration decisions: Case studies from the Barents Sea, Norwegian Continental Shelf, Interpretation , 2( 3), SH55– SH66. Google Scholar CrossRef Search ADS   Gabrielsen P.T., Brevik I., Mittet R., Løseth L.O., 2009. Investigating the exploration potential for 3D CSEM using a calibration survey over the Troll Field, First Break , 27( 6), 67– 75. Grayver A.V., 2015. Parallel three-dimensional magnetotelluric inversion using adaptive finite-element method. Part I: theory and synthetic study, Geophys. J. Int. , 202( 1), 584– 603. Google Scholar CrossRef Search ADS   Grayver A.V., Bürg M., 2014. Robust and scalable 3-D geo-electromagnetic modelling approach using the finite element method, Geophys. J. Int. , 198( 1), 110– 125. Google Scholar CrossRef Search ADS   Grayver A.V., Streich R., Ritter O., 2013. Three-dimensional parallel distributed inversion of CSEM data using a direct forward solver, Geophys. J. Int. , 193( 3), 1432– 1446. Google Scholar CrossRef Search ADS   Gribenko A., Zhdanov M., 2007. Rigorous 3D inversion of marine CSEM data based on the integral equation method, Geophysics , 72( 2), WA73– WA84. Google Scholar CrossRef Search ADS   Holten T., Flekkøy E.G., Singer B., Blixt E.M., Hanssen A., Måløy K.J., 2009. Vertical source, vertical receiver, electromagnetic technique for offshore hydrocarbon exploration, First Break , 27( 5), 89– 93. Hölz S., Swidinsky A., Sommer M., Jegen M., Bialas J., 2015. The use of rotational invariants for the interpretation of marine CSEM data with a case study from the North Alex mud volcano, West Nile Delta, Geophys. J. Int. , 201( 1), 224– 245. Google Scholar CrossRef Search ADS   Jaysaval P., Shantsev D., de la Kethulle de Ryhove S., 2014. Fast multimodel finite-difference controlled-source electromagnetic simulations based on a Schur complement approach, Geophysics , 79( 6), E315– E327. Google Scholar CrossRef Search ADS   Jing C., Green K., Willen D., 2008. CSEM inversion: impact of anisotropy, data coverage, and initial models, in 78th Annual International Meeting , SEG, Expanded Abstracts, pp. 604– 608. Johansen S.E., Amundsen H.E.F., Røsten T., Ellingsrud S., Eidesmo T., Bhuiyan A.H., 2005. Subsurface hydrocarbons detected by electromagnetic sounding, First Break , 23( 3), 31– 36. Key K., 2012. Marine EM inversion using unstructured grids: a 2D parallel adaptive finite element algorithm, in 82nd Annual International Meeting , SEG, Expanded Abstracts, pp. 1– 5. Key K., Du Z., Mattsson J., Mckay A., Midgley J., 2014. Anisotropic 2.5D inversion of towed streamer EM data from three North Sea fields using parallel adaptive finite elements, in 76th Annual International Conference and Exhibition , EAGE, Extended Abstracts, doi:10.3997/2214-4609.20140730. Kilmer M.E., O’Leary D.P., 2001. Choosing regularization parameters in iterative methods for ill-posed problems, SIAM J. Matrix Anal. Appl. , 22( 4), 1204– 1221. Google Scholar CrossRef Search ADS   Kong F.N., Westerdahl H., Ellingsrud S., Eidesmo T., Johansen S., 2002. ‘Seabed Logging’: A possible direct hydrocarbon indicator for deepsea prospects using EM energy, Oil Gas J. , 100, 30– 38. Løseth L.O., Ursin B., Amundsen L., 2007. On the effects of anisotropy in marine CSEM, in 69th Annual International Conference and Exhibition , EAGE, Extended Abstracts, doi:10.3997/2214-4609.201401565. Lu X., Xia C., 2007. Understanding anisotropy in marine CSEM data, in 77th Annual International Meeting , SEG, Expanded Abstracts, pp. 633– 637. Mittet R., Morten J.P., 2013. The marine controlled-source electromagnetic method in shallow water, Geophysics , 78( 2), E67– E77. Google Scholar CrossRef Search ADS   Mohamad S.A., Lorenz L., Hoong L.T., Wei T.K., Chandola S.K., Saadah N., Nazihah F., 2010. A practical example why anisotropy matters—A CSEM case study from South East Asia, in 80th Annual International Meeting , SEG, Expanded Abstracts, pp. 696– 700. Monk P., 2003. Finite Element Methods for Maxwell’s Equations , Oxford University Press, New York. Google Scholar CrossRef Search ADS   Morten J.P., Bjørke A.K., Støren T., 2009. CSEM data uncertainty analysis for 3D inversion, in 79th Annual International Meeting , SEG, Expanded Abstracts, pp. 724– 728. Morten J.P., Roth F., Karlsen S.A., Timko D., Pacurar C., Olsen P.A., Nguyen A.K., Gjengedal J., 2012. Field appraisal and accurate resource estimation from 3D quantitative interpretation of seismic and CSEM data, Leading Edge , 31( 4), 447– 456. Google Scholar CrossRef Search ADS   Newman G.A., Alumbaugh D.L., 1997. Three-dimensional massively parallel electromagnetic inversion —I. Theory, Geophys. J. Int. , 128( 2), 345– 354. Google Scholar CrossRef Search ADS   Newman G.A., Commer M., Carazzone J.J., 2010. Imaging CSEM data in the presence of electrical anisotropy, Geophysics , 75( 2), F51– F61. Google Scholar CrossRef Search ADS   Nocedal J., Wright S.J., 1999. Numerical Optimization , Springer, New York. Google Scholar CrossRef Search ADS   Paige C.C., Saunders M.A., 1982. LSQR: an algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Softw. , 8( 1), 43– 71. Google Scholar CrossRef Search ADS   Pedersen H.T., Hiner M., 2014. Channel play in Foz do Amazonas—exploration and reserve estimate using regional 3D CSEM, First Break , 32( 4), 95– 100. Plessix R.-É., Mulder W.A., 2008. Resistivity imaging with controlled-source electromagnetic data: depth and data weighting, Inverse Probl. , 24( 3), 034012, doi:10.1088/0266-5611/24/3/034012. Google Scholar CrossRef Search ADS   Portniaguine O., Zhdanov M.S., 1999. Focusing geophysical inversion images, Geophysics , 64( 3), 874– 887. Google Scholar CrossRef Search ADS   Ramananjaona C., MacGregor L., 2010. 2.5D inversion of CSEM data in a vertically anisotropic earth, J. Phys. Conf. Ser. , 255( 1), 012004, doi:10.1088/1742-6596/255/1/012004. Google Scholar CrossRef Search ADS   Ramananjaona C., MacGregor L., Andréis D., 2011. Sensitivity and inversion of marine electromagnetic data in a vertically anisotropic stratified earth, Geophys. Prospect. , 59( 2), 341– 360. Google Scholar CrossRef Search ADS   Ren Z., Kalscheuer T., Greenhalgh S., Maurer H., 2013. A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling, Geophys. J. Int. , 194( 2), 700– 718. Google Scholar CrossRef Search ADS   Ren Z., Kalscheuer T., Greenhalgh S., Maurer H., 2014. A finite-element-based domain-decomposition approach for plane wave 3D electromagnetic modeling, Geophysics , 79( 6), E255– E268. Google Scholar CrossRef Search ADS   Rung-Arunwan T., Siripunvaraporn W., 2010. An efficient modified hierarchical domain decomposition for two-dimensional magnetotelluric forward modelling, Geophys. J. Int. , 183( 2), 634– 644. Google Scholar CrossRef Search ADS   Santos R.J., 1996. Equivalence of regularization and truncated iteration for general ill-posed problems, Linear Algebr. Appl. , 236, 25– 33. Google Scholar CrossRef Search ADS   Sasaki Y., 2011. Anisotropic, joint 3D inversion of marine CSEM and MT data, in 81st Annual International Meeting , SEG, Expanded Abstracts, pp. 547– 551. Sasaki Y., 2013. 3D inversion of marine CSEM and MT data: An approach to shallow-water problem, Geophysics , 78( 1), E59– E65. Google Scholar CrossRef Search ADS   Scheunert M., Ullmann A., Afanasjew M., Börner R.-U., Siemon B., Spitzer K., 2016. A cut-&-paste strategy for the 3-D inversion of helicopter-borne electromagnetic data —I. 3-D inversion using the explicit Jacobian and a tensor-based formulation, J. appl. Geophys. , 129, 209– 221. Google Scholar CrossRef Search ADS   Schwalenberg K., Engels M., 2012. Marine controlled source electromagnetic methods for gas hydrate assessment: new instrumentation and first results from the Black Sea test cruise, in Protokoll über das 24. Schmucker-Weidelt-Kolloquium für Elektromagnetische Tiefenforschung , pp. 239– 249, Neustadt a. d. Weinstraße, Germany. Schwarzbach C., 2009. Stability of finite element solutions to Maxwell’s equations in frequency domain, PhD thesis , Technische Universität Bergakademie Freiberg, Germany. Schwarzbach C., Haber E., 2013. Finite element based inversion for time-harmonic electromagnetic problems, Geophys. J. Int. , 193( 2), 615– 634. Google Scholar CrossRef Search ADS   Schwarzbach C., Börner R.-U., Spitzer K., 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics—a marine CSEM example, Geophys. J. Int. , 187( 1), 63– 74. Google Scholar CrossRef Search ADS   Streich R., 2009. 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy, Geophysics , 74( 5), F95– F105. Google Scholar CrossRef Search ADS   Tompkins M.J., 2005. The role of vertical anisotropy in interpreting marine controlled-source electromagnetic data, in 75th Annual International Meeting , SEG, Expanded Abstracts, pp. 514– 517. Weitemeyer K., Gao G., Constable S., Alumbaugh D., 2010. The practical application of 2D inversion to marine controlled-source electromagnetic data, Geophysics , 75( 6), F199– F211. Google Scholar CrossRef Search ADS   Zach J., Bjørke A., Støren T., Maaø F., 2008. 3D inversion of marine CSEM data using a fast finite-difference time-domain forward code and approximate Hessian-based optimization, in 78th Annual International Meeting , SEG, Expanded Abstracts, pp. 614– 618. Zehner B., Börner J.H., Görz I., Spitzer K., 2015. Workflows for generating tetrahedral meshes for finite element simulations on complex geological structures, Comput. Geosci. , 79, 105– 117. Google Scholar CrossRef Search ADS   Zhang Y., Key K., 2016. MARE3DEM: A three-dimensional CSEM inversion based on a parallel adaptive finite element method using unstructured meshes, in 86th Annual International Meeting , SEG, Expanded Abstracts, pp. 1009– 1013. Zhdanov M.S., 2009. Geophysical Electromagnetic Theory and Methods , Elsevier, Amsterdam. Zhdanov M.S., Gribenko A., Čuma M., 2007. Regularized focusing inversion of marine CSEM data using minimum vertical-support stabilizer, in 77th Annual International Meeting , SEG, Expanded Abstracts, pp. 579– 583. Zhdanov M.S., Endo M., Cox L.H., Čuma M., Linfoot J., Anderson C., Black N., Gribenko A.V., 2014a. Three-dimensional inversion of towed streamer electromagnetic data, Geophys. Prospect. , 62( 3), 552– 572. Google Scholar CrossRef Search ADS   Zhdanov M.S., Endo M., Yoon D., Čuma M., Mattsson J., Midgley J., 2014b. Anisotropic 3D inversion of towed-streamer electromagnetic data: Case study from the Troll West Oil Province, Interpretation , 2( 3), SH97– SH113. Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Anisotropic three-dimensional inversion of CSEM data using finite-element techniques on unstructured grids

, Volume 213 (2) – May 1, 2018
17 pages

/lp/ou_press/anisotropic-three-dimensional-inversion-of-csem-data-using-finite-RUSiDhj2bt
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggy029
Publisher site
See Article on Publisher Site

### Abstract

Summary In this paper, we present a recently developed anisotropic 3-D inversion framework for interpreting controlled-source electromagnetic (CSEM) data in the frequency domain. The framework integrates a high-order finite-element forward operator and a Gauss–Newton inversion algorithm. Conductivity constraints are applied using a parameter transformation. We discretize the continuous forward and inverse problems on unstructured grids for a flexible treatment of arbitrarily complex geometries. Moreover, an unstructured mesh is more desirable in comparison to a single rectilinear mesh for multisource problems because local grid refinement will not significantly influence the mesh density outside the region of interest. The non-uniform spatial discretization facilitates parametrization of the inversion domain at a suitable scale. For a rapid simulation of multisource EM data, we opt to use a parallel direct solver. We further accelerate the inversion process by decomposing the entire data set into subsets with respect to frequencies (and transmitters if memory requirement is affordable). The computational tasks associated with each data subset are distributed to different processes and run in parallel. We validate the scheme using a synthetic marine CSEM model with rough bathymetry, and finally, apply it to an industrial-size 3-D data set from the Troll field oil province in the North Sea acquired in 2008 to examine its robustness and practical applicability. Electrical anisotropy, Marine electromagnetics, Inverse theory, Numerical approximations and analysis 1 INTRODUCTION For derisking the exploration of submarine oil and gas resources, the controlled-source electromagnetic (CSEM) method is now regarded as a crucial surveying methodology providing electrical resistivity information on the subsurface (Kong et al.2002; Johansen et al.2005; Constable & Srnka 2007; Holten et al.2009; Anderson & Mattson 2010; Schwalenberg & Engels 2012; Fanavoll et al.2014a,b; Hölz et al.2015). It is complementary to seismics and emerged in the 2000s (Eidesmo et al.2002). A horizontal electric dipole (HED) source can be deployed efficiently and generates a field with feasible polarization to detect thin resistive hydrocarbon-bearing layers. Nowadays, the advances in the development of instruments have greatly improved the data quality and acquisition efficiency such that 3-D dense survey layouts are becoming routine for acquiring large volumes of data (Gabrielsen et al.2009; Pedersen & Hiner 2014). A comprehensive data coverage can enhance the imaging resolution of 3-D structures and thus reduce the ambiguity in data analysis. However, this also creates the need for developing sophisticated interpretation tools such as 3-D inversion algorithms to treat large data sets. An important aspect of interpreting marine CSEM data is to understand the effect of electrical anisotropy (Løseth et al.2007; Ellis et al.2010). The fine layering of marine sediments and asymmetric grain shapes yield a larger resistivity in the direction normal to the deposited layers than parallel to them (Tompkins 2005; Newman et al.2010; Brown et al.2012). In this paper, we assume that the layers have small tilt angles so that the layer normal and parallel directions can be approximated as the vertical and horizontal directions (VTI anisotropic medium, Jing et al.2008). The impact of electrical anisotropy on CSEM data has a complicated dependence on many factors such as water depths and survey geometries. Based on a 1-D scenario with varying degrees of electrical anisotropy, Lu & Xia (2007) illustrated that the inline data are stronger related to the vertical resistivity, while the broadside measurements are more sensitive to the horizontal resistivity. This observation was further elaborated by Ramananjaona et al. (2011) using 1-D forward modeling and inversion. Furthermore, they also observed that inversion of marine CSEM data measured in an anisotropic environment under the assumption of electrical isotropy can lead to a resolved resistivity model with many misleading oscillations (artefacts). Similar conclusions have also been found in both 2-D (Ramananjaona & MacGregor 2010) and 3-D inversion scenarios (Jing et al.2008; Newman et al.2010; Sasaki 2011). Although broadside data are less sensitive to thin resistive layers than inline data, the incorporation of them in an anisotropic inversion is proven to better resolve the electrical anisotropy (Newman et al.2010) and to be useful to reduce the overall data misfit (Mohamad et al.2010). These findings have been promoting the use of full 3-D data sets. Besides allowing for the host medium being electrically anisotropic, Abubakar et al. (2010) and Brown et al. (2012) have discussed the possibility of resolving the electrical anisotropy of a thin reservoir, proposing the addition of other sources such as horizontal magnetic dipoles to enhance the inductive effect for reconstructing the horizontal resistivity. In this paper, however, we only focus on an HED source due to the field data set at hand. In order to interpret marine CSEM data, reported anisotropic 3-D inversion schemes include a parallel non-linear conjugate gradient (CG) optimization approach with a finite-difference forward operator (Newman et al.2010), a quasi-Newton optimizer equipped with a time-domain finite-difference forward code (Zach et al.2008), a Gauss–Newton optimization with a finite-difference forward operator (Sasaki 2013; Amaya et al.2016), a reweighted regularized CG method with an integral equation modeling operator (Zhdanov et al.2014b) and a Schur complement finite-difference forward operator (Jaysaval et al.2014) for serving the anisotropic quasi-Newton optimization approach presented by Zach et al. (2008). In this paper, we introduce a new anisotropic 3-D inversion code which has been implemented using finite-element techniques. Although finite-element methods have received considerable attention for electromagnetic forward modeling (Schwarzbach et al.2011; Ren et al.2013; Ansari & Farquharson 2014; Grayver & Bürg 2014), to our best knowledge, the published inversion results utilizing finite-element approaches were mainly 2-D anisotropic (Key 2012; Key et al.2014) and 3-D isotropic (Schwarzbach & Haber 2013; Grayver 2015; Zhang & Key 2016). Rather than discretizing the forward problem and parametrizing the inverse problem on structured grids, we adopt unstructured mesh discretization. Unstructured meshes offer great flexibility in handling complex geometries to allow for very accurate descriptions of complicated seafloors, reservoir structures, or submetre scale infrastructure like pipelines and casings in the geomodel. Unstructured meshes also make it extremely convenient to carry out local grid refinement. High-order finite-element methods are available for solving the forward problem with high accuracy (Castillo et al.2005; Schwarzbach 2009). We formulate the forward problem by means of the so-called secondary field approach which only simulates scattered fields generated by the discrepancies of electromagnetic parameters from a known background for primary fields (Streich 2009; Schwarzbach et al.2011). Although the total field approach is also implemented in the code, the secondary field approach is favoured here for simulating marine CSEM data. The reason for this is twofold. At first, it is reasonable to assume that the sea surface is flat such that a horizontally stratified model is readily regarded as the background in which field responses can be efficiently obtained. Moreover, the singularity removal technique that the secondary field method applies allows a reduction of the node density around the sources which makes it advantageous over the total field approach when many sources of finite length are considered in the same model (mesh). As marine CSEM data are often processed at only a few frequencies but many transmitter positions (Gabrielsen et al.2009; Weitemeyer et al.2010; Zhdanov et al.2014a), it is desirable to use a parallel direct solver to inexpensively obtain the responses of the transmitters operated at the same frequency. We formulate the inverse problem in an unconstrained least-squares sense, penalized by the smoothness regularization. Although it has been argued that images stabilized by smoothness functionals cannot sufficiently describe the true geological structures (Portniaguine & Zhdanov 1999; Gribenko & Zhdanov 2007; Zhdanov et al.2007), our experiences with this regularization are positive because it is very robust and barely problem-dependent (such as on model parametrization). We solve the discrete inverse problem using an inexact Gauss–Newton approach by virtue of its relatively fast convergence rate. The scheme is implemented in a parallel manner to accelerate the inversion of multisource and multifrequency EM data. Its efficiency and versatility are further considered by providing three alternatives to handle the Jacobian. In the following, we first revisit the finite-element discretization of the forward problem in terms of the secondary electric field. Then, we describe the inverse problem and introduce the inexact Gauss–Newton algorithm. To compute the Jacobian occurring in the Gauss–Newton equation systems, we offer three approaches. The efficiency of each method depends on the survey geometry and problem size and, therefore, the selection is data-driven. Finally, we examine the feasibility of the framework using a synthetic marine example with rough bathymetry and the Troll 3-D CSEM field data set acquired in 2008 in the North Sea (Gabrielsen et al.2009). 2 FORWARD PROBLEM We assume a harmonic time dependence e−iωt to simulate EM fields in the space $$\mathbf {r}\in \mathbb {R}^3$$ induced by an electric source with a current density Js = Js(r, ω). We construct a bounded computational domain $$\Omega \in \mathbb {R}^3$$ with its outer boundary $$\partial \Omega \in \mathbb {R}^3$$ far enough from the source to properly characterize the true field behaviour. The associated boundary value problem consists of the curl–curl equation of the electric field and a simple Dirichlet boundary condition imposed at the outer boundary ∂Ω to give a unique solution   \begin{eqnarray} \nabla \times \,({{{\mu }}}^{-1}\nabla \times \,{\mathbf {E}}) - i \omega \boldsymbol{\tilde{\varsigma }}{\mathbf {E}}= i \omega \mathbf {J}_s \quad {\rm in} \quad \Omega , \end{eqnarray} (1a)  \begin{eqnarray} \mathbf {n}\times {\mathbf {E}}= \mathbf {0} \quad {\rm on} \quad \partial \Omega , \end{eqnarray} (1b) where n denotes the unit normal vector to the outer boundary, μ is the magnetic permeability which is assumed to be isotropic and allowed to vary throughout the model here and $$\boldsymbol{\tilde{\varsigma }}$$ denotes the tensor complex electrical conductivity. In a medium with transverse isotropy about a vertical axis, $$\boldsymbol{\tilde{\varsigma }}$$ is characterized by a 3 × 3 diagonal matrix   $$\boldsymbol{\tilde{\varsigma }}= \left(\begin{array}{ccc}\sigma _h - i \omega \epsilon _h &0 &0 \\ 0 & \sigma _h - i \omega \epsilon _h &0 \\ 0 & 0 &\sigma _v - i \omega \epsilon _v \end{array} \right),$$ (2)where σh and σv mean the horizontal and vertical conductivities, respectively, and εh and εv denote the horizontal and vertical electrical permittivities, respectively. With the working frequencies in the range between about 0.1 and 10 Hz (Constable 2010) and by assigning a small conductivity value such as 10−7 Sm−1 to the air layer (Jaysaval et al.2014), the marine CSEM method works in the quasi-static limit where σ ≫ ωε holds. However, we preserve the full-wave expression to generalize the applicability of our code. To avoid the numerical difficulty in directly simulating the rapidly decaying field around the transmitter, the total electric field E = Ep + Es is split into a known primary field Ep due to the impressed current Js in a background model ($${{\mu }}_0, \boldsymbol{\tilde{\varsigma }}_0$$) and a secondary electric field Es to be solved for. In doing so, the mesh node density around the source can be relaxed compared to the one for formulation (1). Moving the terms including Ep to the right-hand side of eq. (1) attains the boundary value problem in terms of Es (Schwarzbach et al.2011)   \begin{eqnarray} \nabla \!\times ({{{\mu }}}^{-1}\nabla \!\times {\mathbf {E}^s}) \!-\! i \omega \boldsymbol{\tilde{\varsigma }}{\mathbf {E}^s}\!=\! \nabla \!\times (\delta {{{\mu }}}^{-1}\nabla \!\times {\mathbf {E}^p}) \!-\! i \omega \delta \boldsymbol{\tilde{\varsigma }}{\mathbf {E}^p}\ \ {\rm in } \quad \Omega ,\nonumber\\ \end{eqnarray} (3a)  \begin{eqnarray} \mathbf {n} \times {\mathbf {E}^s}= \mathbf {0} \ \ {\rm on }\quad \partial \Omega , \end{eqnarray} (3b) where the right-hand side of eq. (3a) now acts as the secondary sources and $$\delta {{{\mu }}}^{-1}= {{{\mu }}}^{-1}_0 - {{{\mu }}}^{-1}$$ and $$\delta \boldsymbol{\tilde{\varsigma }}= \boldsymbol{\tilde{\varsigma }}_0 - \boldsymbol{\tilde{\varsigma }}$$ are the parameter discrepancies between the background and total models. We solve boundary value problem (3) using the Nédélec curl-conforming finite-element method on unstructured grids (Castillo et al.2005; Schwarzbach et al.2011). The computational domain Ω is decomposed into a set of non-overlapping tetrahedra first. The unknown function Es is then approximated in the curl-conforming finite-element function space which is spanned by small-support basis functions   $$\mathbf {E}^s(\mathbf {r}) \approx \sum _{i=1}^{N_{\mathrm{dof}}} E^s_i \mathbf {N}_i(\mathbf {r}),$$ (4)where Ndof denotes the total number of degrees of freedom, Ni means the ith Nédélec curl-conforming basis function (Monk 2003), and $$E^s_i$$ denotes the ith discrete unknown coefficient in the mesh. After assembling all the finite-element matrices, solving the original boundary value problem (3) is reduced to solving the following system of linear equations   $$\mathbf {K} \mathbf {x} = \mathbf {F},$$ (5)where $$\mathbf {x} = (E_1^s, \ldots , E^s_{N_{\mathrm{dof}}})^T \in \mathbb {C}^{N_{\mathrm{dof}}}$$ denotes the solution vector, $$\mathbf {K} \in \mathbb {C}^{N_{\mathrm{dof}}\times N_{\mathrm{dof}}}$$ is the system matrix which is very sparse and symmetric, and $$\mathbf {F} \in \mathbb {C}^{N_{\mathrm{dof}}}$$ represents the secondary source vector which possesses non-zero entries where the electrical or magnetic parameters of the model differ from the background. In detail, the elements of K and F read (Schwarzbach et al.2011)   \begin{eqnarray} {K}_{ij} \!=\! \int _{\Omega } \nabla \times \,\mathbf {N}_i \cdot ({{{\mu }}}^{-1}\nabla \times \,\mathbf {N}_j) {d\Omega }- i \omega \int _{\Omega } \mathbf {N}_i \cdot \left(\boldsymbol{\tilde{\varsigma }}\mathbf {N}_j \right) {d\Omega }, \end{eqnarray} (6a)and   \begin{eqnarray} F_{i} \!=\! \int _{\Omega } \nabla \times \,\mathbf {N}_i \cdot (\delta {{{\mu }}}^{-1}\nabla \times {\mathbf {E}^p}) {d\Omega }\!-\! i \omega \int _{\Omega } \mathbf {N}_i \cdot \left(\delta \boldsymbol{\tilde{\varsigma }}{\mathbf {E}^p}\right) {d\Omega }. \end{eqnarray} (6b) Due to the sparsity of the matrix K, the system of eq. (5) can be solved either iteratively or directly. Considering that realistic marine surveys often yield thousands of transmitter positions, we prefer a parallel direct solver which might outperform an iterative approach for moderate-size problems because the solutions for the transmitters at one frequency can be computed rapidly by merely forward and backward substitutions once the system matrix is decomposed and the factors are stored. 3 INVERSE PROBLEM We discretize the forward and inverse problems on the same unstructured mesh which allows for incorporating all transmitters, frequencies, and receivers, and, at the same time, provides a reasonable inversion resolution. Assuming that the inversion domain Ωm coincides with or is a part of the forward model Ω: Ωm ⊆ Ω and includes Nm non-overlapping subdomains, we parametrize the conductivity model using the pixel-based method (e.g. Commer & Newman 2009). In each subdomain, the horizontal (σh) and vertical (σv) conductivities are decoupled and treated as independent unknown parameters   $$\left( \begin{array}{c}\sigma _h(\mathbf {r})\\ \sigma _v(\mathbf {r}) \end{array} \right) \approx \sum _{i=1}^{N_m} \psi _i(\mathbf {r}) \left( \begin{array}{c}\sigma _{h_i}\\ \sigma _{v_i} \end{array}\right),$$ (7)where $$\sigma _{h_i}$$ and $$\sigma _{v_i}$$ refer to the unknown horizontal and vertical conductivity coefficients (degrees of freedom) with respect to the ith inversion subdomain $$\Omega _{m_i}$$, respectively, and ψi(r) denotes the corresponding basis function which is non-trivial and equals unity only over $$\Omega _{m_i}\!$$. To ensure that the reconstructed electrical conductivity is always positive and finite, we enforce the bounded model transformation with a specified lower and upper conductivity bound set to σa and σb, respectively (Commer & Newman 2008). Instead of directly recovering the electrical conductivity, the transformed parameter m will be solved for first. The true conductivity is then obtained by the inverse model transformation. The forward and inverse transformations are given by   \begin{eqnarray} \left( \begin{array}{c}m_{h_i}\\ m_{v_i} \end{array} \right) \!=\!\left( \begin{array}{c}\log (\sigma _{h_i} \!-\! \sigma _a)\!-\! \log (\sigma _b \!-\! \sigma _{h_i})\\ \log (\sigma _{v_i} \!-\! \sigma _a)\!-\! \log (\sigma _b \!-\! \sigma _{v_i}) \end{array} \right); \sigma _a<\sigma _{h_i}, \sigma _{v_i} < \sigma _b,\nonumber\\ \end{eqnarray} (8a)  \begin{eqnarray} \left( \begin{array}{c}\sigma _{h_i}\\ \sigma _{v_i} \end{array}\right) = \left( \begin{array}{c}\frac{\sigma _a+\sigma _be^{m_{h_i}}}{1+e^{m_{h_i}}}\\ \frac{\sigma _a+\sigma _be^{m_{v_i}}}{1+e^{m_{v_i}}} \end{array}\right); -\infty <m_{h_i}, m_{v_i} < \infty . \end{eqnarray} (8b) It should be mentioned that the horizontal and vertical conductivity bounds (namely σa and σb in the above equations) do not necessarily have to be identical. As most typically in the VTI assumption, the horizontal conductivity is higher than the vertical one, a practical constraint for the horizontal conductivity may have a strict lower bound and a relaxed upper bound. In contrast, a loose lower bound and a tight upper bound can be more effective for the vertical conductivity. However, in the absence of any a priori information, we have used the same bounds for all the examples in the paper. Given the parametrized model and assuming that the measurements $$\mathbf {d}^\mathrm{obs}\in \mathbb {C}^{N_d}$$ are collected at Nf frequencies, Ns transmitters (at $${\mathbf {r}^s_{i}}, i = 1,\ldots ,N_s$$), and Nr receivers (at $${\mathbf {r}^r_{i}}, i = 1,\ldots , N_r$$) for Nt field components (labeled in a vector $$\boldsymbol{\Theta }$$), namely Nd = Nf × Ns × Nr × Nt, we seek a parameter model m which minimizes the following regularized discrete objective functional   \begin{eqnarray} \Phi (\mathbf {m}) &=& \Phi _d(\mathbf {m}) + \lambda _h \Phi _{m_h}(\mathbf {m})+\lambda _v \Phi _{m_v}(\mathbf {m}) \nonumber\\ &=& \frac{1}{2} \Vert \mathbf {W}(\mathbf {d}^\mathrm{obs}- \mathbf {d}^{\mathrm{syn}}(\mathbf {m}^{}))\Vert ^2 + \frac{\lambda _h}{2}\Vert \mathbf {L}_h(\mathbf {m}^{}_h-\mathbf {m}^{\mathrm{ref}}_h)\Vert ^2 \nonumber\\ &&+ \frac{\lambda _v}{2}\Vert \mathbf {L}_v(\mathbf {m}^{}_v-\mathbf {m}^{\mathrm{ref}}_v)\Vert ^2, \end{eqnarray} (9)where, in the data misfit functional, that is, the first term on the right-hand side, dsyn(m) represents the predicted data vector and W is a diagonal matrix storing the data weights. The determination of proper data weights is crucial for a successful inversion of CSEM data because the long-offset fields are orders of magnitude weaker than the short-offset ones. To treat the short- and long-offset data somehow equally, a commonly used data weight involves the standard deviation of a datum (denoted by δ|dobs| for the amplitude and δϕobs for the phase throughout the rest of the paper) (e.g. Grayver et al.2013; Schwarzbach & Haber 2013). We construct the data weighting matrix for the synthetic data set according to Schwarzbach & Haber (2013) which incorporates the consequences of the number of data, the relative error in the data amplitude, a scaling factor for different data components and the noise floor. As for the field data, a different weighting is employed to account for the rotation uncertainties which will be illustrated in Section 4.2. In the VTI inversion, the model vector $$\mathbf {m} = \left({\mathbf {m}^{}_h}^T, {\mathbf {m}^{}_v}^T\right)^T$$ is composed of the horizontal model parameter vector $$\mathbf {m}^{}_h = ({m}^{}_{h_{1}}, \ldots , {m}^{}_{h_{N_m}})^T \in \mathbb {R}^{N_m}$$ and the vertical model parameter vector $$\mathbf {m}^{}_v = ({m}^{}_{v_{1}}, \ldots , {m}^{}_{v_{N_m}})^T \in \mathbb {R}^{N_m}$$ which results in twice as many model parameters as the isotropic inversion; λh and λv in eq. (9) are the trade-off or regularization parameters; $$\Phi _{m_h}(\mathbf {m})$$ and $$\Phi _{m_v}(\mathbf {m})$$ are the model misfit functionals for the horizontal and vertical parameter models, respectively, which are chosen to be the discrete form of the smoothness stabilizer measuring the gradient of the model discrepancy in the l2 norm   $$\Phi _m(m) = \frac{1}{2}\int _{\Omega _m} \vert \nabla (m-m^{\mathrm{ref}}) \vert ^2 {d\Omega }\qquad {\rm in} \ {\Omega _m}.$$ (10)On an unstructured mesh, the discrete form of the above functional can be obtained by resorting to the lowest order divergence-conforming Nédélec finite-element method (Schwarzbach & Haber 2013). Lh and Lv in eq. (9) are identical representing the resulting smoothness matrix which is independent of the model parameters; mrefh and mrefv are the reference models for $$\mathbf {m}^{}_h$$ and $$\mathbf {m}^{}_v$$, respectively. To handle the data topology, a global indexing matrix $$\boldsymbol{\mathcal {I}} \in \mathbb {R}^{{N_d}\times 4}$$ is built. Its row number refers to the global index of a data point (i). In each row, there are four elements, in turn representing the global index of the frequency ($${\mathcal {I}_{i1}} \in [1,N_f]$$), transmitter ($${\mathcal {I}_{i2}} \in [1,N_s]$$), receiver ($${\mathcal {I}_{i3}} \in [1, N_r]$$) and data type ($${\mathcal {I}_{i4}}$$). To calculate the ith data point $$d^\mathrm{syn}_i(\mathbf {m}^{})$$ on the current parameter model m, the real physical model $${\boldsymbol{\sigma }^{}}$$ must be obtained via the inverse model transformation (8b) first. Second, the discrete forward problem (5) with the model $${\boldsymbol{\sigma }^{}}$$ due to the frequency-transmitter pair $$(\omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})$$ is solved for the solution $${{\mathbf {x}}}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})$$. Finally, $$d^\mathrm{syn}_i(\mathbf {m}^{})$$ can be obtained by combing the primary field computed analytically in the background model with the interpolated finite-element solution   \begin{eqnarray} d^\mathrm{syn}_i(\mathbf {m}^{}) &=& d^p_i(\omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}}, {\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, \Theta _{{\mathcal {I}_{i4}}}) \nonumber\\ &&+ \mathbf {q}^T(\omega _{\mathcal {I}_{i1}}\!,{\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, \Theta _{{\mathcal {I}_{i4}}}) {{\mathbf {x}}}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}}), \end{eqnarray} (11)where $$\mathbf {q}(\omega _{\mathcal {I}_{i1}}\!,{\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, \Theta_{{\mathcal {I}_{i4}}}) \in \mathbb {C}^{N_{\mathrm{dof}}}$$ is the (complex) interpolation vector built to evaluate the field component $$\Theta_{{\mathcal {I}_{i4}}}$$ at the frequency $$\omega _{\mathcal {I}_{i1}}$$ and receiver $${\mathbf {r}^r_{{\mathcal {I}_{i3}}}}$$. It does not depend on the transmitter position and is constructed according to the interpolation rules of the curl-conforming finite elements and thus includes only geometric weights and frequency-related values for specific data types. For instance, if the datum is the x-component of the electric field, the specific expression of q reads   $$\mathbf {q}(\omega _{\mathcal {I}_{i1}}\!,{\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, E_x) = \left\lbrace N_{1,x}({\mathbf {r}^r_{{\mathcal {I}_{i3}}}}),\ldots , N_{\mathrm{dof},x}({\mathbf {r}^r_{{\mathcal {I}_{i3}}}})\right\rbrace ^T.$$ (12)Because the basis functions have a small support, q turns out to be very sparse. 3.1 An inexact Gauss–Newton framework To find a minimizer of the objective functional (9), we employ an inexact Gauss–Newton approach. In terms of the data Jacobian, the gradient ∇Φ(m) and the Gauss–Newton Hessian HGN(m) of the functional can be expressed as   \begin{eqnarray} {\nabla \Phi (\mathbf {m}) = }\nonumber\\ &&\, && \Re \left\lbrace \begin{array}{c}-[\mathbf {W}{\mathbf {J}_h(\mathbf {m}^{})}]^H \mathbf {W}[\mathbf {d}^\mathrm{obs}- \mathbf {d}^{\mathrm{syn}}(\mathbf {m}^{})] + \lambda _h \mathbf {L}_h^T\mathbf {L}_h(\mathbf {m}^{}_h - \mathbf {m}^{\mathrm{ref}}_h) \\ \quad -[\mathbf {W}{\mathbf {J}_v(\mathbf {m}^{})}]^H \mathbf {W}[\mathbf {d}^\mathrm{obs}- \mathbf {d}^{\mathrm{syn}}(\mathbf {m}^{})] + \lambda _v \mathbf {L}_v^T\mathbf {L}_v(\mathbf {m}^{}_v - \mathbf {m}^{\mathrm{ref}}_v) \end{array} \right\rbrace , \nonumber\\ \end{eqnarray} (13a)  \begin{eqnarray} {{\mathbf {H}^{\rm GN}(\mathbf {m})}}\nonumber\\ && {{ \!=\! \Re \!\left\lbrace \begin{array}{c}\left[\mathbf {W}{\mathbf {J}_h(\mathbf {m}^{})}\right]^H [\mathbf {W}{\mathbf {J}_h(\mathbf {m}^{})}] \!+\! \lambda _h \mathbf {L}_h^T\mathbf {L}_h \quad\ [\mathbf {W}{\mathbf {J}_h(\mathbf {m}^{})}]^H [\mathbf {W}{\mathbf {J}_v(\mathbf {m}^{})}] \\ \quad \left[\mathbf {W}{\mathbf {J}_v(\mathbf {m}^{})}\right]^H [\mathbf {W}{\mathbf {J}_h(\mathbf {m}^{})}] \quad\ [\mathbf {W}{\mathbf {J}_v(\mathbf {m}^{})}]^H [\mathbf {W}{\mathbf {J}_v(\mathbf {m}^{})}]\!+\!\lambda _v \mathbf {L}_v^T\mathbf {L}_v \end{array} \right\rbrace ,}}\nonumber\\ \end{eqnarray} (13b) where $${\mathbf {J}_h(\mathbf {m}^{})} \in \mathbb {C}^{N_d\times N_m}$$ and $${\mathbf {J}_v(\mathbf {m}^{})} \in \mathbb {C}^{N_d \times N_m}$$ are the data Jacobian components with respect to $$\mathbf {m}^{}_h$$ and $$\mathbf {m}^{}_v$$, respectively,   \begin{eqnarray} \left. {\mathbf {J}_h(\mathbf {m}^{})} \right|_{ij} &=& \frac{\partial d^\mathrm{syn}_i(\mathbf {m}^{})}{\partial {m}^{}_{h_{j}}}=-\mathbf {q}^T(\omega _{\mathcal {I}_{i1}}\!,{\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, \Theta _{{\mathcal {I}_{i4}}}) \mathbf {K}^{-1}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}) \nonumber \\ \quad \times\left[\frac{\partial \mathbf {K}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}})}{\partial {\sigma }^{}_{h_{j}}} {{\mathbf {x}}}({\boldsymbol{\sigma }},\omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})\right.\nonumber\\ \quad \left. - \frac{\partial {\mathbf {F}^{{}}}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})}{\partial {\sigma }^{}_{h_{j}}}\right]{\frac{\partial {\sigma }^{}_{h_{j}}}{\partial {m}^{}_{h_{j}}}}, \end{eqnarray} (14a)  \begin{eqnarray} \left. {\mathbf {J}_v(\mathbf {m}^{})} \right|_{ij} &=& \frac{\partial d^\mathrm{syn}_i(\mathbf {m}^{})}{\partial {m}^{}_{v_{j}}} = -\mathbf {q}^T(\omega _{\mathcal {I}_{i1}}\!,{\mathbf {r}^r_{{\mathcal {I}_{i3}}}}\!, \Theta _{{\mathcal {I}_{i4}}}) \mathbf {K}^{-1}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}) \nonumber \\ \quad \times \left[\frac{\partial \mathbf {K}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}})}{\partial {\sigma }^{}_{v_{j}}} {{\mathbf {x}}}({\boldsymbol{\sigma }},\omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})\right.\nonumber\\ \quad -\left. \frac{\partial {\mathbf {F}^{{}}}({\boldsymbol{\sigma }}, \omega _{\mathcal {I}_{i1}}\!, {\mathbf {r}^s_{{\mathcal {I}_{i2}}}})}{\partial {\sigma }^{}_{v_{j}}}\right]{\frac{\partial {\sigma }^{}_{v_{j}}}{\partial {m}^{}_{v_{j}}}}. \end{eqnarray} (14b) We note that the mixed partial derivatives of the stabilizing functionals vanish because the regularizations on $$\mathbf {m}^{}_h$$ and $$\mathbf {m}^{}_v$$ are decoupled. Moreover, we see that the off-diagonal matrices of the Hessian (13b) actually couple the solutions through anisotropy. The Gauss–Newton method minimizes the non-linear functional (9) through solving the following linearized subproblem iteratively   $$\mathbf {H}^{\rm GN}(\mathbf {m}^{k}) \left( \begin{array}{c}\delta \mathbf {m}^{k}_h\\ \delta \mathbf {m}^{k}_v \end{array} \right) = -\nabla \Phi (\mathbf {m}^{k}),$$ (15)where k denotes the kth Gauss–Newton iteration. For practical 3-D CSEM inverse problems, this dense equation system has to be solved iteratively which only requires matrix-vector operations. To distinguish the approach for solving eq. (15) from the one for eq. (5), we denote the former by the inner solver. The normal eq. (15) might become very ill-conditioned for large-scale inverse problems with a non-uniform parametrization. Rather than resorting to the conventional CG method, a numerically more stable strategy is to solve the system (15) in the linear least-squares sense using some least-squares solver such as LSQR (Paige & Saunders 1982) as it only operates on J instead of JHJ  \begin{eqnarray} {\left( \begin{array}{c@\quad c}\Re \lbrace \mathbf {W}{\mathbf {J}_h(\mathbf {m}^{k})}\rbrace & \Re \lbrace \mathbf {W}{\mathbf {J}_v(\mathbf {m}^{k})}\rbrace \\ \Im \lbrace \mathbf {W}{\mathbf {J}_h(\mathbf {m}^{k})}\rbrace & \Im \lbrace \mathbf {W}{\mathbf {J}_v(\mathbf {m}^{k})}\rbrace \\ \sqrt{\lambda _h^k} \mathbf {L}_{h} & \mathbf {0} \\ \mathbf {0} & \sqrt{\lambda _v^k} \mathbf {L}_{v}\\ \end{array} \right) \left( \begin{array}{c}\delta \mathbf {m}^{k}_h\\ \delta \mathbf {m}^{k}_v \end{array} \right)} \nonumber\\ =\left( \begin{array}{cc}\Re \lbrace \mathbf {W}[\mathbf {d}^{\mathrm{obs}}- \mathbf {d}^{\mathrm{syn}}(\mathbf {m}^{k})]\rbrace \\ \Im \lbrace \mathbf {W}[\mathbf {d}^{\mathrm{obs}}- \mathbf {d}^{\mathrm{syn}}(\mathbf {m}^{k})]\rbrace \\ -\sqrt{\lambda _h^k} \mathbf {L}_{h}(\mathbf {m}^{k}_h-\mathbf {m}^{\mathrm{ref}}_h) \\ -\sqrt{\lambda _v^k} \mathbf {L}_{v}(\mathbf {m}^{k}_v-\mathbf {m}^{\mathrm{ref}}_v) \end{array}\right). \end{eqnarray} (16) For large-scale problems, on the one hand, solving eq. (16) exactly is impractical. On the other hand, it is rather unnecessary to obtain an exact solution especially when the current model mk is far away from a (local) minimum. Some studies suggest that regularization is implicitly applied when the solving procedure is terminated at early iterations (Santos 1996; Kilmer & O’Leary 2001; Grayver et al.2013). The truncation of the iteration is problem-dependent and many authors have found some useful stopping criteria by empirical investigations to explain the observations to the level of measurement errors (e.g. Grayver et al.2013; Schwarzbach & Haber 2013). The Gauss–Newton method with an early termination when solving the linearized subproblem (16) is termed the inexact Gauss–Newton method (Nocedal & Wright 1999). To prevent the Gauss–Newton method from being divergent, an inexact line search algorithm according to the Armijo condition is augmented to steer the decrease in the objective functional (Nocedal & Wright 1999). 3.2 The data Jacobian In a Gauss–Newton-based inversion scheme, preparing the full data Jacobian at each iteration is the most time-consuming ingredient besides solving the forward problem. Determining an efficient strategy for treating the Jacobian is problem-dependent. To accommodate most of the cases, we have provided three alternatives in this framework to deal with the Jacobian. For small- or moderate-size problems, the Jacobian can be constructed explicitly. According to eq. (14), this will lead to a memory cost of $$\mathcal {O}(N_d\times N_m\times 2)$$ by carrying out Nf × Nr × Nt additional forward runs. With an explicit Jacobian, operating the Jacobian-vector products becomes very efficient which allows the inner solver to carry out more iterations for a more accurate solution (Scheunert et al.2016). In contrast to the explicit storage, a widely used strategy is to treat the Jacobian implicitly which works for an iterative solution of the Gauss–Newton equation system (15). If the inner iterative process is terminated after Nitr iterations, 2 × Nf × Ns × Nitr additional forward runs are needed to obtain a model update δmk. This approach is memory-efficient, but time-consuming because the required number of forward simulations increases linearly with the number of inner iterations Nitr. Besides the above two approaches, we introduce a third method which we call the semi-explicit strategy. It acts as a compromise between the explicit and implicit methods because only one part of the Jacobian is stored explicitly which requires forward runs. At each Gauss–Newton iteration, we first compute and store the vectors qTK − 1 occurring in the Jacobian (i.e. the solution of the adjoint problem at each receiver, see eq. 14) by exploiting the symmetry of the system matrix K, which needs a memory of $$\mathcal {O}(N_f \times N_r \times N_t \times N_{\mathrm{dof}})$$ and, like the explicit storage, Nf × Nr × Nt additional forward runs (Newman & Alumbaugh 1997). Any Jacobian-vector multiplication Jp can be done by first carrying out a sparse matrix-vector product $$(\frac{\partial \mathbf {K}}{\partial \boldsymbol{\sigma }} \mathbf {x} - \frac{\partial \mathbf {F}}{\partial \boldsymbol{\sigma }})\frac{\partial \boldsymbol{\sigma }}{\partial \mathbf {m}} \mathbf {p}$$ and then multiplying the pre-computed row vectors by the resulting vector. The inverse procedure applies when operating the transpose of the data Jacobian. The memory cost in the semi-explicit method might be less demanding than that in the explicit storage if there are only tens to hundreds of receiver positions which can be the case in marine CSEM applications (Gabrielsen et al.2009). The efficiency of Jacobian-vector operations compared to that in the explicit method largely relies on the relation between Nm × 2 and Ndof. For high-order curl-conforming basis functions, the latter can be much larger than the former one, while they might be comparable when the lowest order basis function is considered. Nonetheless, a notable advantage of this strategy is that the memory requirement is not affected by the number of model parameters. This makes it competitive for large-scale anisotropic 3-D CSEM inverse problems in which the number of model parameters increases linearly with the dimension of anisotropy. 4 TEST AND APPLICATION In this section, we first examine the new anisotropic inversion scheme using a synthetic marine model with rough bathymetry. Based on this, we apply it to interpret the Troll field data acquired in 2008 by Statoil ASA and EMGS ASA in the North Sea (Gabrielsen et al.2009). All computations were conducted on a moderately parallel machine equipped with 4 Intel(R) Xeon(R) 2.2 GHz CPUs with 8 cores per CPU and 512 GB of RAM. To accelerate the inversion process, we have implemented a parallel scheme typically designed for marine CSEM data at a few frequencies and many transmitter positions. First, the computational model is built and discretized considering all frequencies, transmitters and receivers. Then, the entire data set is decomposed with respect to frequencies. As solving the system of linear equations for the forward problem with a large number of transmitters can be a significant burden when using a direct solver, each data subset can be further decomposed with respect to transmitters to accelerate the solution phase. Each computational task which consists of computing the Jacobian-vector products, the gradient and misfit of the objective function for the corresponding data subset is then run in parallel using the Message Passing Interface specification. Inside each process, the finite-element equation systems are solved using MUMPS (Amestoy et al.2001). For a direct solver, the decomposition with respect to transmitters is optional when sufficient memory is available. For the synthetic example, we decompose the data set with respect to frequencies. For the field data set, the parallel decomposition is also performed with respect to the transmitters. The inner solver LSQR is terminated when either the maximum number of iterations (120 for the synthetic example and 60 for the field experiment) is achieved or the tolerance 10−2 is satisfied. The regularization parameters λh and λv are chosen to be equal unless otherwise stated. We terminate the inversion process once the data misfit does not decrease significantly. 4.1 A synthetic marine CSEM model The model consists of an air layer of 10−7 Sm−1, a sea water column of 3.7037 Sm−1, an undulating seafloor, a subsurface with a horizontal conductivity of 1 Sm−1 and a vertical conductivity of 0.25 Sm−1, and two resistive anomalies of 0.01 Sm−1. Considering a right-handed Cartesian coordinate system with the z-axis pointing downwards, the flat sea surface is situated at z = 0 m. The depth of the sea water in the survey area varies from about z = 352 m to about z = 916 m (see the depth contour in Fig. 1a). The smaller resistive body has a dimension of 1500 m × 2000 m × 100 m along the x-, y- and z-directions with the centre being located at (2750 m, −500 m, 1400 m). The larger body has a size of 3000 m × 3000 m × 100 m and is centred at (− 1500 m, 0 m, 1500 m). Figure 1. View largeDownload slide Survey geometry over the undulating seafloor. (a) It displays the depth contour of the seafloor topography and the dashed black lines outline the horizontal extents of the anomalies. Receivers are denoted by black triangles and transmitters by red dots. (b) It depicts the vertical cross-section of the subsurface along transmitter line Tx03 showing the depth extents of the two anomalous bodies. (c) It highlights the bathymetry variation along each towline Tx01 to Tx04. Figure 1. View largeDownload slide Survey geometry over the undulating seafloor. (a) It displays the depth contour of the seafloor topography and the dashed black lines outline the horizontal extents of the anomalies. Receivers are denoted by black triangles and transmitters by red dots. (b) It depicts the vertical cross-section of the subsurface along transmitter line Tx03 showing the depth extents of the two anomalous bodies. (c) It highlights the bathymetry variation along each towline Tx01 to Tx04. We consider a dense CSEM survey including four transmitter lines named Tx01 to Tx04. As shown in Fig. 1, 44 seabed receivers (black triangles) are arranged in an 11 × 4 rectangular array with an x-directed spacing of 900 m and a y-directed spacing of 1000 m. Each receiver is equipped with two orthogonal electric antennae for recording the horizontal components of the electric field $$E_{x^{\prime }}$$ and $$E_{y^{\prime }}$$. Due to the bathymetry, the sensors generally have a non-zero pitch (right-handed rotation around the y-axis, positive bow up) which varies from −15.8° to 12.0° as shown in Fig. 2. Consequently, $$E_{x^{\prime }}$$ and $$E_{y^{\prime }}$$ might not be horizontally oriented. We do not rotate the components into the global coordinate system but denote and compute them in each receiver’s local coordinates using the exact rotation angles. The vertical position of the receivers varies from around z = 438 m to around z = 794 m. An x-directed HED of 200 m length is towed about 30 m above each x-directed receiver line ranging from x = −8100 to 8100 m with an acquisition separation of 300 m (red dots in Fig. 1). This gives rise to a total number of 220 transmitter positions. The transmission frequencies were 0.25 and 0.75 Hz. Both inline and offline data within the transmitter–receiver range between 500 and 9000 m were recorded yielding 33 792 complex measurements in total. Figure 2. View largeDownload slide Receiver pitch. Figure 2. View largeDownload slide Receiver pitch. Since the number of transmitters (220) is larger than the product of the number of receivers (44) and the number of data types (2), we apply the reciprocity principle (Zhdanov 2009). By doing this, the real receiver antennae will become the computational transmitters. Correspondingly, the computational receivers are now placed where the real transmitters are located to record the component along the transmitter’s orientation. This leads to 44 x΄-directed and 44 y΄-directed computational electric transmitters (both in the local coordinate system) and 220 x-directed electric receivers with a length of 200 m. We first prepare the model for generating the synthetic data set. It contains the two anomalies embedded in the rugged subsurface with a size of about 60000 m × 60000 m × 15000 m, the sea water column and an air layer with a thickness of 70000 m on top. The whole computational domain is decomposed into 746662 tetrahedra. The second-order finite-element method is employed to solve the forward problem yielding 4724736 degrees of freedom. The background model is chosen to be the air–water half-space (the same has been used for the inversion). Although not shown here, we have checked the modeling accuracy on this mesh by comparing to the simulation result of a planar layer. The synthetic data were normalized by the dipole moment (200 Am) and 3 per cent Gaussian noise was added. To avoid the inverse crime, we set up a new mesh for the inversion. The size of the computational domain stays the same but the geometries of the two anomalies are excluded (Fig. 3). Since the success of a pixel-based parametrization strategy greatly relies on a proper discretization of the inversion domain, we refine the subsurface below the receiver grid uniformly from about z = 1100 to 1700 m to ensure a sufficient resolution of the thin resistors. The computational domain is discretized into 672 154 tetrahedra. The forward operator is based on the second-order finite-element discretization giving 4255 588 degrees of freedom. We consider the entire subsurface which includes 427 509 cells as the inversion domain. We set the lower and upper conductivity bounds to 0.002 and 10 Sm−1, respectively, in the model transformation (8). For this setup, the explicit storage of the anisotropic Jacobian requires 33 792 × 427509 × 2 × 16 B ≈ 430.5 GB of memory which is not affordable on the current computing system. In contrast, the semi-explicit storage strategy merely needs 2 × 220 × 4255588 × 16 B ≈ 27.9 GB of memory. Table 1 lists a comparison of all three methods for treating the Jacobian in terms of the run times and memory cost. From this example, it can be seen that the semi-explicit storage outperforms the implicit storage in terms of speed and the explicit storage in terms of memory requirement. Figure 3. View largeDownload slide Model discretization for synthetic data generation (upper panel) and the inversion (lower panel) at (a) and (d) x = 0 m, (b) and (e) y = 0 m and (c) and (f) z = 1400 m. Figure 3. View largeDownload slide Model discretization for synthetic data generation (upper panel) and the inversion (lower panel) at (a) and (d) x = 0 m, (b) and (e) y = 0 m and (c) and (f) z = 1400 m. Table 1. The performance of different strategies for treating the Jacobian.   Explicit  Implicit  Semi-explicit  One Jacobian-vector product (s)  –  170.5  13.5  Storage (GB)  430.5   0   27.9     Explicit  Implicit  Semi-explicit  One Jacobian-vector product (s)  –  170.5  13.5  Storage (GB)  430.5   0   27.9   View Large Moreover, it is worthwhile to mention that, in the case of the semi-explicit storage, the memory cost for storing the forward solution vectors and the vectors qTK − 1 is the same as used without applying the reciprocity principle. However, the use accelerates the computation of the primary fields and the assembly of the source terms of the forward problem (88 electric dipoles with reciprocity in comparison with 220 electric line sources without reciprocity). The first inversion experiment (E0) starts with a homogeneous isotropic subsurface of 0.5 Sm−1 which differs from the true anisotropic background (σh = 1 Sm−1 and σv = 0.25 Sm−1). We choose the reference model to be the same as the starting model (also for the following four tests). The initial value of the regularization parameters $$\lambda _h^0$$ and $$\lambda _v^0$$ is 0.000164. We decrease them by a factor of 0.6 at each iteration. Fig. 4 shows the slices of the final inversion results. For the sake of illustration, all inversion results are displayed in terms of the electrical resistivity. We observe from this example that a comprehensive data coverage with both inline and azimuthal data can determine the electrical anisotropy in the sediment. The maximum value of the resolved vertical resistivity in the reservoirs is about 30 Ωm compared to the true value 100 Ωm. But the reconstructed resistors are also thicker as seen in the vertical sections. Therefore, it results in a reasonable anomalous transverse resistance (about 11495 Ωm2 for the large body and 4389 Ωm2 for the small one) compared to the true value (9600 Ωm2 for both targets). Anomalous transverse resistance is the quantity that determines the CSEM response. It is also seen from Fig. 4 that the model was updated in the regions with weak sensitivity which is the effect of the smoothness operator. Figs 5 and 6 display the relative amplitude and phase residuals in terms of the data uncertainty. The comparison with the ‘exact’ relative data residuals for the noise-free synthetic data [denoted by (d*, ϕ*)] indicates that the predicted model fits the observed data well in the order of magnitude of the noise. Figure 4. View largeDownload slide Slices of the resolved horizontal (left-hand panel) and vertical (right-hand panel) resistivity models in the experiment E0 at (a) and (b) z = 1480 m, (c) and (d) z = 1360 m, along towline (e) and (f) Tx01, (g) and (h) Tx02, (i) and (j) Tx03 and (k) and (l) Tx04. Figure 4. View largeDownload slide Slices of the resolved horizontal (left-hand panel) and vertical (right-hand panel) resistivity models in the experiment E0 at (a) and (b) z = 1480 m, (c) and (d) z = 1360 m, along towline (e) and (f) Tx01, (g) and (h) Tx02, (i) and (j) Tx03 and (k) and (l) Tx04. Figure 5. View largeDownload slide Relative amplitude misfit in the experiment E0 at (a) the 0th iteration and (b) the 13th iteration in comparison with (c) the true misfit for noise-free data (d*, ϕ*). Transmitters and receivers are indexed by firstly walking along the positive y-direction and then the positive x-direction. Figure 5. View largeDownload slide Relative amplitude misfit in the experiment E0 at (a) the 0th iteration and (b) the 13th iteration in comparison with (c) the true misfit for noise-free data (d*, ϕ*). Transmitters and receivers are indexed by firstly walking along the positive y-direction and then the positive x-direction. Figure 6. View largeDownload slide Relative phase misfit in the experiment E0 at (a)  the 0th iteration and (b)  the 13th iteration in comparison with (c) the true misfit for noise-free data (d*, ϕ*). Figure 6. View largeDownload slide Relative phase misfit in the experiment E0 at (a)  the 0th iteration and (b)  the 13th iteration in comparison with (c) the true misfit for noise-free data (d*, ϕ*). Figure 7. View largeDownload slide The development of the (a) data misfit function and (b)  stabilizing functions for different inversion conditions in the synthetic model studies. Figure 7. View largeDownload slide The development of the (a) data misfit function and (b)  stabilizing functions for different inversion conditions in the synthetic model studies. To further investigate the dependence of the inversion results on the starting model, the number of inner iterations Nitr, and the model constraints, we carry out another four experiments. To examine whether the current inversion result is affected by insufficiently solving each Gauss–Newton system of equations (16) or a loose model constraint, we perform the first additional model study (E1) by setting Nitr = 150 and a stricter model bound of σ ∈ (0.005, 2) Sm−1. The second (E2) and third (E3) additional experiments examine the influence of different starting models. E2 considers a homogeneous subsurface of σh = σv = 1 Sm−1 as the starting model with Nitr = 120 and a bound constraint σ ∈ (0.002, 10) Sm−1. E3 utilizes the same inversion conditions for E2 except for the true background model as the initial guess. The fourth additional experiment (E4) is based on experiment E3 but using Nitr = 300. Figs 7a and 7b plot the development of the data misfit function Φd and regularization functions $$\Phi _{m_h}$$ and $$\Phi _{m_v}$$ for the above five experiments, respectively. Despite the use of different inversion conditions, the five tests have converged to the same level: E0 converges to a total data misfit of 0.5520, E1 is 0.5531, E2 is 0.5493, E3 is 0.5515 and E4 is 0.5511. Fig. 8 depicts the resistivities along vertical pseudo-boreholes as a function of depth at the horizontal receiver locations of (−1800 m, 0 m) and (2700 m, −1000 m). We see that a large number of LSQR iterations in E4 does not significantly improve the inversion result indicating Nitr = 120 is sufficient for the current example. The stronger model constraint used in E1 seems to enhance the resolved anomalies compared to other setups. Figure 8. View largeDownload slide Resistivities along a vertical borehole through the anisotropic inversion models at the horizontal positions of (a) and (b) (2700 m, −1000 m) and (c) and (d)( − 1800 m, 0 m). Figure 8. View largeDownload slide Resistivities along a vertical borehole through the anisotropic inversion models at the horizontal positions of (a) and (b) (2700 m, −1000 m) and (c) and (d)( − 1800 m, 0 m). We also note that the horizontal resistivity of the two anomalies (100 Ωm) is hardly resolved with this survey configuration. This phenomenon is caused by the guided wave excited in the thin resistive layer. Due to the continuity equation, a very strong vertical electric field is generated inside the resistor. The magnitude of the horizontal electric field is in comparison negligible which prevents the horizontal resistivity from being resolved (Ramananjaona et al.2011; Mittet & Morten 2013). Moreover, Brown et al. (2012) pointed out that using an HED source, the horizontal resistivity of a thin resistive reservoir cannot be reconstructed once its value is greater than that of the background. This finding corresponds well with the observation here. 4.2 Troll field The Troll field is located in a mature area in the North Sea where early CSEM surveys showed that the CSEM method can be a powerful exploration technique due to the prominent variation in EM fields caused by the large-scale resistive gas reservoir. Some authors have reported the imaging results of the CSEM data set acquired in 2003 using a single receiver line across the Troll West Gas Province (TWGP, e.g. Commer & Newman 2008; Plessix & Mulder 2008). Here, we utilize a more comprehensive CSEM data set collected in 2008 using a dense 3-D acquisition grid by Statoil ASA and EMGS ASA to map the resistivity distribution of this area (Gabrielsen et al.2009; Morten et al.2012; Jaysaval et al.2014). The whole survey contains six north–south towlines being around 1.25 km apart, three east–west towlines with a large line spacing of about 3.8 km, and 53 stationary seabed detectors aiming to focus on the relatively small-scale Troll West Oil Province (TWOP) while covering only a small part of the TWGP (see Gabrielsen et al.2009). Both inline and azimuthal data were recorded simultaneously and processed at frequencies 0.25, 0.75 and 1.25 Hz. The inline measurements from the westernmost receiver line were excluded because of the large impact of pipelines. For a clear illustration, we only show the processed layout in Fig. 9 by excluding the westernmost receiver line. The processed data at hand are the rotated horizontal components of the electric field Ex and Ey with Ex being along the corresponding towline direction and Ey being perpendicular to that. The data uncertainties of the rotated field components are related to the data uncertainties of the non-rotated field components through error propagation analysis and the data weights are constructed based on that according to Morten et al. (2009). Figure 9. View largeDownload slide The adapted Troll field survey layout in 2008 (Gabrielsen et al.2009). The colour map describes the variation of the water depth. Black rectangles indicate seabed receivers and red lines represent towlines. The TWOP is outlined in white and TWGP in black. Figure 9. View largeDownload slide The adapted Troll field survey layout in 2008 (Gabrielsen et al.2009). The colour map describes the variation of the water depth. Black rectangles indicate seabed receivers and red lines represent towlines. The TWOP is outlined in white and TWGP in black. For the inversion, we select the inline and azimuthal data from the five north–south and three east–west towlines at frequencies 0.25 and 0.75 Hz within the transmitter–receiver range between 1200 and 9000 m. The data subset includes 157 748 data from 44 Ex and 44 Ey receivers and 38 349 transmitter positions. Again, we apply the reciprocity principle resulting in 88 computational seabed transmitters and 38 349 towline-directed electric receivers of 270.5 m length. The bathymetry over the survey area (Fig. 9) is rather smooth and small compared to the survey scale. The water depth varies from about 300 m in the southeast to about 360 m in the northeast. Although previous studies have assumed that the bathymetry over this region does not influence the measurements much (Newman et al.2010), we incorporated it into the model to highlight the flexibility of unstructured meshes and the finite-element approach. We build and visualize the inversion models using a coordinate system with x along the north direction, y along the east direction and z pointing downwards. The sea surface is at z = 0 m. The subsurface with a topographic seafloor has a size of about 70000 m × 66000 m × 30000 m, overlain by sea water with a conductivity of 3.7037 Sm−1, and an air layer of 10−7 Sm−1. We prepared distinct meshes for isotropic and anisotropic inversions, respectively. This is reasonable because our preliminary study showed that many small-scale structures were generated close to the seafloor in the isotropic inversion. Therefore, it shall require a denser discretization in that case to rule out the possibility that this is caused by the inaccuracy of the forward simulation. In order to ensure a reasonable resolution for both oil and gas fields, we refine the subdomain which laterally covers the horizontal extents of the oil and gas fields and vertically extends from the sea bottom to about z = 2700 m. The vertical spacing of the cells in the refined region is around 100 m. Outside the subdomain, the mesh gets coarser towards the outer boundaries. The regions around the transmitters and receivers are locally refined with a minimum edge length of around 0.4 m to ensure sufficient accuracy. For anisotropic inversions, the whole computational domain is decomposed into 4882 262 tetrahedra, 2269 297 of which are distributed in the sea water layer and 2337 595 in the subsurface. The cells close to the seafloor are further refined in the isotropic case. The resulting denser mesh consists of 5152 631 cells, 2401 811 of which are in the sea water and 2463 768 in the subsurface. For all experiments, we employ the first-order finite-element method to solve the forward problem with an air–water half-space as the background model. The whole subsurface is considered as the inversion domain. The model constraint is set to σ ∈ (0.002, 5) Sm−1. Other inversion setups and statistics are listed in Table 2. For this large-scale configuration, neither the explicit (too many model parameters) nor the semi-explicit (too many computational receivers) storage of the data Jacobian is feasible on our current computing system. Therefore, the Jacobian is treated implicitly. Each inversion task is decomposed into 2 × 2 (number of partitions over frequencies × number of partitions over transmitters) inversion jobs. With a direct solver, the drawback of this decomposition is the repeated storage of the factor of the system matrix at the same frequency. However, if the memory requirement is permitted, a scalability factor of (nearly) 2 is observed compared to a 2 × 1 decomposition. Table 2. Inversion conditions for the three Troll experiments.   Isotropic  Anisotropic (inline)  Anisotropic (inline and azimuthal)  Mesh size  5152 631  4882 262  4882 262  Nm  2463 768  2337 595  2337 595  Nd  16 088  16 088  157 748  Ndof  5985 264  5672 052  5672 052  $$\lambda _h^0$$ and $$\lambda _v^0$$  0.000252  0.000207  0.0000255  mrefh and mrefv  σ = 0.5 Sm−1  σh = σv = 0.5 Sm−1  σh = σv = 0.5 Sm−1  Run time (hour)  59.1  48.2  67.4  Matrix factor (GB)  67.3  62.6  62.6    Isotropic  Anisotropic (inline)  Anisotropic (inline and azimuthal)  Mesh size  5152 631  4882 262  4882 262  Nm  2463 768  2337 595  2337 595  Nd  16 088  16 088  157 748  Ndof  5985 264  5672 052  5672 052  $$\lambda _h^0$$ and $$\lambda _v^0$$  0.000252  0.000207  0.0000255  mrefh and mrefv  σ = 0.5 Sm−1  σh = σv = 0.5 Sm−1  σh = σv = 0.5 Sm−1  Run time (hour)  59.1  48.2  67.4  Matrix factor (GB)  67.3  62.6  62.6  View Large First, we carry out an isotropic inversion of the inline data from the selected data subset. The left-hand panel of Fig. 10 shows the slices of the resolved resistivity model at the depth of z = 1500 m and along the five north–south towlines. Although both oil and gas fields are indicated, there exist many small-scale anomalies which likely are imaging artefacts if compared with a priori geostructures. This poses the question of whether electrical anisotropy is present. To demonstrate this assumption, secondly, we run an anisotropic inversion of the same inline data subset and display the recovered vertical resistivity model in the right-hand column of Fig. 10. The comparison shows that the anisotropic inversion renders much fewer artefacts around the receivers and near the seafloor than the isotropic inversion. Moreover, the anisotropic resistivity model becomes geologically more plausible and the horizontal extents of the resistivity anomalies coincide well with the outlined a priori regions. The oscillating layers observed in the overburden in the isotropic inversion, which can to some extent represent an effect like anisotropy, have largely disappeared from the anisotropic inversion. To show an example how the models fit the measurements, the relative data residuals at the common offset of 2500 m are checked for both inversions. They are positioned at the common midpoint of the corresponding transmitters and receivers. It is seen from Figs 11 and 12 that the anisotropic resistivity model yields a better data fit than the isotropic resistivity model. Based on the above observations, electrical anisotropy seems to be significant in this data set. Figure 10. View largeDownload slide Slices of the isotropic (left-hand panel) and anisotropic (right-hand panel) Troll resistivity models characterized by merely inline data at (a) and (b)  z = 1500 m, along towline (c) and (d)  Tx002, (e) and (f) Tx003, (g) and (h) Tx004, (i) and (j) Tx005 and (k) and (l) Tx006. The solid black and white lines in (a) and (b) outline the horizontal extents of the Troll gas and oil fields, respectively, and black triangles indicate receiver positions. Figure 10. View largeDownload slide Slices of the isotropic (left-hand panel) and anisotropic (right-hand panel) Troll resistivity models characterized by merely inline data at (a) and (b)  z = 1500 m, along towline (c) and (d)  Tx002, (e) and (f) Tx003, (g) and (h) Tx004, (i) and (j) Tx005 and (k) and (l) Tx006. The solid black and white lines in (a) and (b) outline the horizontal extents of the Troll gas and oil fields, respectively, and black triangles indicate receiver positions. Figure 11. View largeDownload slide Relative amplitude misfit from the isotropic (left-hand panel) and anisotropic (right-hand panel) inversions of the Troll inline data in the common-midpoint domain at the offset of 2500 m. Figure 11. View largeDownload slide Relative amplitude misfit from the isotropic (left-hand panel) and anisotropic (right-hand panel) inversions of the Troll inline data in the common-midpoint domain at the offset of 2500 m. Figure 12. View largeDownload slide Relative phase misfit from the isotropic (left-hand panel) and anisotropic (right-hand panel) inversions of the Troll inline data in the common-midpoint domain at the offset of 2500 m. Figure 12. View largeDownload slide Relative phase misfit from the isotropic (left-hand panel) and anisotropic (right-hand panel) inversions of the Troll inline data in the common-midpoint domain at the offset of 2500 m. Finally, we run another anisotropic inversion by additionally incorporating the azimuthal data. The slices of the imaged horizontal and vertical resistivity models are displayed in Fig. 13 and the selected data residuals in Fig. 14. It is seen through comparing Figs 13 and 10 that, despite the fact that the anisotropic inversion of purely inline data has yielded the main characteristics of the Troll oil and gas reservoirs, the incorporation of the azimuthal data has noticeably enhanced the reconstruction of both fields: (i) both the TWOP and TWGP are now indicated more clearly; (ii) their lateral extents are enlarged and (iii) the up-tilted part of the gas reservoir, which is strongly observed in the isotropic inversion while slightly shown in the corresponding anisotropic case, is now pushed down. According to the geological description and the inversion results shown by Gabrielsen et al. (2009) and Morten et al. (2012), the oil and gas fields have been resolved at the correct depth. The depth extents of both oil and gas fields are a bit too large ranging from about z = 1250 to about 1650 m for the oil field and from about z = 1100 to about 1800 m for the gas field. This is probably caused by the limited resolution of the data, the use of the smoothness regularization, and the relatively large vertical spacing of the cells (about 100 m). Furthermore, as shown in Fig. 15, the two anisotropic inversion processes steadily converge to a smaller overall data misfit than the isotropic inversion procedure. Figure 13. View largeDownload slide Slices of the anisotropic Troll resistivity model characterized by both inline and azimuthal data at depth (a) and (b) z = 1500 m, along towline (c) and (d) Tx002, (e) and (f) Tx003, (g) and (h) Tx004, (i) and (j) Tx005 and (k) and (l) Tx006. Figure 13. View largeDownload slide Slices of the anisotropic Troll resistivity model characterized by both inline and azimuthal data at depth (a) and (b) z = 1500 m, along towline (c) and (d) Tx002, (e) and (f) Tx003, (g) and (h) Tx004, (i) and (j) Tx005 and (k) and (l) Tx006. Figure 14. View largeDownload slide Relative amplitude (left-hand panel) and phase (right-hand panel) misfits from the anisotropic inversion of the full Troll data subset in the common-midpoint domain at the offset of 2500 m. Figure 14. View largeDownload slide Relative amplitude (left-hand panel) and phase (right-hand panel) misfits from the anisotropic inversion of the full Troll data subset in the common-midpoint domain at the offset of 2500 m. Figure 15. View largeDownload slide The development of the Troll (a) data and (b)  regularization misfits at each Gauss–Newton iteration. Note that the curve for $$\Phi _{m_h}$$ is absent in the isotropic case. Figure 15. View largeDownload slide The development of the Troll (a) data and (b)  regularization misfits at each Gauss–Newton iteration. Note that the curve for $$\Phi _{m_h}$$ is absent in the isotropic case. 5 CONCLUSIONS In this paper, we have presented an anisotropic 3-D inversion framework based on a finite-element discretization on unstructured grids for interpreting CSEM data in the frequency domain. The complexity of finite-element modeling compared to that of finite-difference modeling is balanced by its flexibility in incorporating different types of meshes. The freedom that comes with an unstructured mesh allows to describe very complicated geometries and avoids interpolation problems that are typical for rectilinear mesh approaches. This advantage becomes even more attractive when simulating large 3-D CSEM data sets with many transmitter positions involved in complex geometries. However, the application of finite-element methods is, to some extent, impeded by the effort of generating a quality mesh for a very complex model. To deal with this issue, some promising workflows have been reported in which complicated 3-D surfaces are represented in terms of non-uniform rational basis splines (Börner et al.2015; Zehner et al.2015). The resulting free surface can therefore be flexibly triangulated to avoid mesh degeneration caused by the description of a 3-D undulating surface through interpolating a fixed triangulation. The inversion scheme has incorporated two different formulations, namely by using the total and secondary field approaches for the forward problem, respectively, and offered three different strategies for treating the data Jacobian for a general application. Moreover, we have achieved a parallelism that allows us to tackle an industry-scale 3-D data set on a moderate-scale parallel machine. We restricted ourselves to marine CSEM applications in this paper. We picked the secondary field approach and validated the functionalities of the inversion scheme using a close-to-real-world synthetic marine CSEM scenario and the realistic Troll field data set from a 3-D dense survey grid in 2008. In this framework, we chose to use the parallel direct solver MUMPS for solving the system of linear equations for the multisource forward problem. Although direct solvers are oftentimes criticized because of high memory demand and poor scalability for large-scale problems, the Troll field example with about 5.7 million degrees of freedom has shown that the performance of the direct solver is still competitive. The memory requirement in a problem of such a scale is not only determined by the direct solver. The memory cost for storing the solutions, the sources (the right-hand side of the forward equation system), and some other auxiliary vectors for Jacobian-related operations is also non-trivial. However, when the problem size gets even larger, one might need to further optimize the forward solver by considering methods such as domain decomposition strategies (Rung-Arunwan & Siripunvaraporn 2010; Bakr et al.2013; Ren et al.2014) or/and fast convergent iterative approaches to relieve the memory request for storing the factor of the entire system matrix. We put this into our future research task list. Both the synthetic and real model studies have shown that the horizontal resistivity of thin resistive reservoirs cannot be resolved. The lack of sensitivity to the horizontal resistivity might be overcome by adding measurements due to other sources which generate a TE-mode dominant pattern in thin resistive anomalies. The significance of taking electrical anisotropy into account for marine CSEM data has been addressed by comparing the isotropic and anisotropic inversion results of the Troll field data. Acknowledgements We are grateful to EMGS ASA for providing the Troll field data set and for the permission to publish the results. FW would like to thank Christoph Schwarzbach for offering the finite-element forward operator. The work was partially funded by the China Scholarship Council (no. 2009637116). And finally, we thank editor Gary Egbert, the two anonymous reviewers, and Colin Farquharson for their insightful comments. REFERENCES Abubakar A., Liu J., Li M., Habashy T.M., MacLennan K., 2010. Sensitivity study of multi-sources receivers CSEM data for TI-anisotropy medium using 2.5D forward and inversion algorithm, in 72nd Annual International Conference and Exhibition , EAGE, Extended Abstracts, doi:10.3997/2214-4609.201400664 Amaya M., Morten J.P., Boman L., 2016. A low-rank approximation for large-scale 3D controlled-source electromagnetic Gauss-Newton inversion, Geophysics , 81( 3), E211– E225. Google Scholar CrossRef Search ADS   Amestoy P.R., Duff I.S., L’Excellent J.-Y., Koster J., 2001. A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM J. Matrix Anal. Appl. , 23( 1), 15– 41. Google Scholar CrossRef Search ADS   Anderson C., Mattson J., 2010. An integrated approach to marine electromagnetic surveying using a towed streamer and source, First Break , 28( 5), 71– 75. Ansari S., Farquharson C.G., 2014. 3D finite-element forward modeling of electromagnetic data using vector and scalar potentials and unstructured grids, Geophysics , 79( 4), E149– E165. Google Scholar CrossRef Search ADS   Bakr S.A., Pardo D., Mannseth T., 2013. Domain decomposition Fourier finite element method for the simulation of 3D marine CSEM measurements, J. Comput. Phys. , 255, 456– 470. Google Scholar CrossRef Search ADS   Börner J.H., Wang F., Weißflog J., Bär M., Görz I., Spitzer K., 2015. Multi-method virtual electromagnetic experiments for developing suitable monitoring designs: A fictitious CO2 sequestration scenario in Northern Germany, Geophys. Prospect. , 63( 6), 1430– 1449. Google Scholar CrossRef Search ADS   Brown V., Hoversten M., Key K., Chen J., 2012. Resolution of reservoir scale electrical anisotropy from marine CSEM data, Geophysics , 77( 2), E147– E158. Google Scholar CrossRef Search ADS   Castillo P., Rieben R., White D., 2005. FEMSTER: an object-oriented class library of high-order discrete differential forms, ACM Trans. Math. Softw. , 31( 4), 425– 457. Google Scholar CrossRef Search ADS   Commer M., Newman G.A., 2008. New advances in three-dimensional controlled-source electromagnetic inversion, Geophys. J. Int. , 172( 2), 513– 535. Google Scholar CrossRef Search ADS   Commer M., Newman G.A., 2009. Three-dimensional controlled-source electromagnetic and magnetotelluric joint inversion, Geophys. J. Int. , 178( 3), 1305– 1316. Google Scholar CrossRef Search ADS   Constable S., 2010. Ten years of marine CSEM for hydrocarbon exploration, Geophysics , 75( 5), 75A67– 75A81. Google Scholar CrossRef Search ADS   Constable S., Srnka L.J., 2007. An introduction to marine controlled-source electromagnetic methods for hydrocarbon exploration, Geophysics , 72( 2), WA3– WA12. Google Scholar CrossRef Search ADS   Eidesmo T., Ellingsrud S., MacGregor L.M., Constable S., Sinha M.C., Johansen S.E., Kong F.N., Westerdahl H., 2002. Sea Bed Logging (SBL), a new method for remote and direct identification of hydrocarbon filled layers in deepwater areas, First Break , 20( 3), 144– 152. Ellis M.H., Sinha M., Parr R., 2010. Role of fine-scale layering and grain alignment in the electrical anisotropy of marine sediments, First Break , 28( 9), 49– 57. Fanavoll S., Gabrielsen P.T., Ellingsrud S., 2014a. The impact of CSEM on exploration decisions and seismic: two case studies from the Barents Sea, First Break , 32( 11), 105– 110. Fanavoll S., Gabrielsen P.T., Ellingsrud S., 2014b. CSEM as a tool for better exploration decisions: Case studies from the Barents Sea, Norwegian Continental Shelf, Interpretation , 2( 3), SH55– SH66. Google Scholar CrossRef Search ADS   Gabrielsen P.T., Brevik I., Mittet R., Løseth L.O., 2009. Investigating the exploration potential for 3D CSEM using a calibration survey over the Troll Field, First Break , 27( 6), 67– 75. Grayver A.V., 2015. Parallel three-dimensional magnetotelluric inversion using adaptive finite-element method. Part I: theory and synthetic study, Geophys. J. Int. , 202( 1), 584– 603. Google Scholar CrossRef Search ADS   Grayver A.V., Bürg M., 2014. Robust and scalable 3-D geo-electromagnetic modelling approach using the finite element method, Geophys. J. Int. , 198( 1), 110– 125. Google Scholar CrossRef Search ADS   Grayver A.V., Streich R., Ritter O., 2013. Three-dimensional parallel distributed inversion of CSEM data using a direct forward solver, Geophys. J. Int. , 193( 3), 1432– 1446. Google Scholar CrossRef Search ADS   Gribenko A., Zhdanov M., 2007. Rigorous 3D inversion of marine CSEM data based on the integral equation method, Geophysics , 72( 2), WA73– WA84. Google Scholar CrossRef Search ADS   Holten T., Flekkøy E.G., Singer B., Blixt E.M., Hanssen A., Måløy K.J., 2009. Vertical source, vertical receiver, electromagnetic technique for offshore hydrocarbon exploration, First Break , 27( 5), 89– 93. Hölz S., Swidinsky A., Sommer M., Jegen M., Bialas J., 2015. The use of rotational invariants for the interpretation of marine CSEM data with a case study from the North Alex mud volcano, West Nile Delta, Geophys. J. Int. , 201( 1), 224– 245. Google Scholar CrossRef Search ADS   Jaysaval P., Shantsev D., de la Kethulle de Ryhove S., 2014. Fast multimodel finite-difference controlled-source electromagnetic simulations based on a Schur complement approach, Geophysics , 79( 6), E315– E327. Google Scholar CrossRef Search ADS   Jing C., Green K., Willen D., 2008. CSEM inversion: impact of anisotropy, data coverage, and initial models, in 78th Annual International Meeting , SEG, Expanded Abstracts, pp. 604– 608. Johansen S.E., Amundsen H.E.F., Røsten T., Ellingsrud S., Eidesmo T., Bhuiyan A.H., 2005. Subsurface hydrocarbons detected by electromagnetic sounding, First Break , 23( 3), 31– 36. Key K., 2012. Marine EM inversion using unstructured grids: a 2D parallel adaptive finite element algorithm, in 82nd Annual International Meeting , SEG, Expanded Abstracts, pp. 1– 5. Key K., Du Z., Mattsson J., Mckay A., Midgley J., 2014. Anisotropic 2.5D inversion of towed streamer EM data from three North Sea fields using parallel adaptive finite elements, in 76th Annual International Conference and Exhibition , EAGE, Extended Abstracts, doi:10.3997/2214-4609.20140730. Kilmer M.E., O’Leary D.P., 2001. Choosing regularization parameters in iterative methods for ill-posed problems, SIAM J. Matrix Anal. Appl. , 22( 4), 1204– 1221. Google Scholar CrossRef Search ADS   Kong F.N., Westerdahl H., Ellingsrud S., Eidesmo T., Johansen S., 2002. ‘Seabed Logging’: A possible direct hydrocarbon indicator for deepsea prospects using EM energy, Oil Gas J. , 100, 30– 38. Løseth L.O., Ursin B., Amundsen L., 2007. On the effects of anisotropy in marine CSEM, in 69th Annual International Conference and Exhibition , EAGE, Extended Abstracts, doi:10.3997/2214-4609.201401565. Lu X., Xia C., 2007. Understanding anisotropy in marine CSEM data, in 77th Annual International Meeting , SEG, Expanded Abstracts, pp. 633– 637. Mittet R., Morten J.P., 2013. The marine controlled-source electromagnetic method in shallow water, Geophysics , 78( 2), E67– E77. Google Scholar CrossRef Search ADS   Mohamad S.A., Lorenz L., Hoong L.T., Wei T.K., Chandola S.K., Saadah N., Nazihah F., 2010. A practical example why anisotropy matters—A CSEM case study from South East Asia, in 80th Annual International Meeting , SEG, Expanded Abstracts, pp. 696– 700. Monk P., 2003. Finite Element Methods for Maxwell’s Equations , Oxford University Press, New York. Google Scholar CrossRef Search ADS   Morten J.P., Bjørke A.K., Støren T., 2009. CSEM data uncertainty analysis for 3D inversion, in 79th Annual International Meeting , SEG, Expanded Abstracts, pp. 724– 728. Morten J.P., Roth F., Karlsen S.A., Timko D., Pacurar C., Olsen P.A., Nguyen A.K., Gjengedal J., 2012. Field appraisal and accurate resource estimation from 3D quantitative interpretation of seismic and CSEM data, Leading Edge , 31( 4), 447– 456. Google Scholar CrossRef Search ADS   Newman G.A., Alumbaugh D.L., 1997. Three-dimensional massively parallel electromagnetic inversion —I. Theory, Geophys. J. Int. , 128( 2), 345– 354. Google Scholar CrossRef Search ADS   Newman G.A., Commer M., Carazzone J.J., 2010. Imaging CSEM data in the presence of electrical anisotropy, Geophysics , 75( 2), F51– F61. Google Scholar CrossRef Search ADS   Nocedal J., Wright S.J., 1999. Numerical Optimization , Springer, New York. Google Scholar CrossRef Search ADS   Paige C.C., Saunders M.A., 1982. LSQR: an algorithm for sparse linear equations and sparse least squares, ACM Trans. Math. Softw. , 8( 1), 43– 71. Google Scholar CrossRef Search ADS   Pedersen H.T., Hiner M., 2014. Channel play in Foz do Amazonas—exploration and reserve estimate using regional 3D CSEM, First Break , 32( 4), 95– 100. Plessix R.-É., Mulder W.A., 2008. Resistivity imaging with controlled-source electromagnetic data: depth and data weighting, Inverse Probl. , 24( 3), 034012, doi:10.1088/0266-5611/24/3/034012. Google Scholar CrossRef Search ADS   Portniaguine O., Zhdanov M.S., 1999. Focusing geophysical inversion images, Geophysics , 64( 3), 874– 887. Google Scholar CrossRef Search ADS   Ramananjaona C., MacGregor L., 2010. 2.5D inversion of CSEM data in a vertically anisotropic earth, J. Phys. Conf. Ser. , 255( 1), 012004, doi:10.1088/1742-6596/255/1/012004. Google Scholar CrossRef Search ADS   Ramananjaona C., MacGregor L., Andréis D., 2011. Sensitivity and inversion of marine electromagnetic data in a vertically anisotropic stratified earth, Geophys. Prospect. , 59( 2), 341– 360. Google Scholar CrossRef Search ADS   Ren Z., Kalscheuer T., Greenhalgh S., Maurer H., 2013. A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling, Geophys. J. Int. , 194( 2), 700– 718. Google Scholar CrossRef Search ADS   Ren Z., Kalscheuer T., Greenhalgh S., Maurer H., 2014. A finite-element-based domain-decomposition approach for plane wave 3D electromagnetic modeling, Geophysics , 79( 6), E255– E268. Google Scholar CrossRef Search ADS   Rung-Arunwan T., Siripunvaraporn W., 2010. An efficient modified hierarchical domain decomposition for two-dimensional magnetotelluric forward modelling, Geophys. J. Int. , 183( 2), 634– 644. Google Scholar CrossRef Search ADS   Santos R.J., 1996. Equivalence of regularization and truncated iteration for general ill-posed problems, Linear Algebr. Appl. , 236, 25– 33. Google Scholar CrossRef Search ADS   Sasaki Y., 2011. Anisotropic, joint 3D inversion of marine CSEM and MT data, in 81st Annual International Meeting , SEG, Expanded Abstracts, pp. 547– 551. Sasaki Y., 2013. 3D inversion of marine CSEM and MT data: An approach to shallow-water problem, Geophysics , 78( 1), E59– E65. Google Scholar CrossRef Search ADS   Scheunert M., Ullmann A., Afanasjew M., Börner R.-U., Siemon B., Spitzer K., 2016. A cut-&-paste strategy for the 3-D inversion of helicopter-borne electromagnetic data —I. 3-D inversion using the explicit Jacobian and a tensor-based formulation, J. appl. Geophys. , 129, 209– 221. Google Scholar CrossRef Search ADS   Schwalenberg K., Engels M., 2012. Marine controlled source electromagnetic methods for gas hydrate assessment: new instrumentation and first results from the Black Sea test cruise, in Protokoll über das 24. Schmucker-Weidelt-Kolloquium für Elektromagnetische Tiefenforschung , pp. 239– 249, Neustadt a. d. Weinstraße, Germany. Schwarzbach C., 2009. Stability of finite element solutions to Maxwell’s equations in frequency domain, PhD thesis , Technische Universität Bergakademie Freiberg, Germany. Schwarzbach C., Haber E., 2013. Finite element based inversion for time-harmonic electromagnetic problems, Geophys. J. Int. , 193( 2), 615– 634. Google Scholar CrossRef Search ADS   Schwarzbach C., Börner R.-U., Spitzer K., 2011. Three-dimensional adaptive higher order finite element simulation for geo-electromagnetics—a marine CSEM example, Geophys. J. Int. , 187( 1), 63– 74. Google Scholar CrossRef Search ADS   Streich R., 2009. 3D finite-difference frequency-domain modeling of controlled-source electromagnetic data: Direct solution and optimization for high accuracy, Geophysics , 74( 5), F95– F105. Google Scholar CrossRef Search ADS   Tompkins M.J., 2005. The role of vertical anisotropy in interpreting marine controlled-source electromagnetic data, in 75th Annual International Meeting , SEG, Expanded Abstracts, pp. 514– 517. Weitemeyer K., Gao G., Constable S., Alumbaugh D., 2010. The practical application of 2D inversion to marine controlled-source electromagnetic data, Geophysics , 75( 6), F199– F211. Google Scholar CrossRef Search ADS   Zach J., Bjørke A., Støren T., Maaø F., 2008. 3D inversion of marine CSEM data using a fast finite-difference time-domain forward code and approximate Hessian-based optimization, in 78th Annual International Meeting , SEG, Expanded Abstracts, pp. 614– 618. Zehner B., Börner J.H., Görz I., Spitzer K., 2015. Workflows for generating tetrahedral meshes for finite element simulations on complex geological structures, Comput. Geosci. , 79, 105– 117. Google Scholar CrossRef Search ADS   Zhang Y., Key K., 2016. MARE3DEM: A three-dimensional CSEM inversion based on a parallel adaptive finite element method using unstructured meshes, in 86th Annual International Meeting , SEG, Expanded Abstracts, pp. 1009– 1013. Zhdanov M.S., 2009. Geophysical Electromagnetic Theory and Methods , Elsevier, Amsterdam. Zhdanov M.S., Gribenko A., Čuma M., 2007. Regularized focusing inversion of marine CSEM data using minimum vertical-support stabilizer, in 77th Annual International Meeting , SEG, Expanded Abstracts, pp. 579– 583. Zhdanov M.S., Endo M., Cox L.H., Čuma M., Linfoot J., Anderson C., Black N., Gribenko A.V., 2014a. Three-dimensional inversion of towed streamer electromagnetic data, Geophys. Prospect. , 62( 3), 552– 572. Google Scholar CrossRef Search ADS   Zhdanov M.S., Endo M., Yoon D., Čuma M., Mattsson J., Midgley J., 2014b. Anisotropic 3D inversion of towed-streamer electromagnetic data: Case study from the Troll West Oil Province, Interpretation , 2( 3), SH97– SH113. Google Scholar CrossRef Search ADS   © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations