TY - JOUR
AU - Kuvshinov,, Alexey
AB - Summary 3-D interpretation of electromagnetic (EM) data of different origin and scale becomes a common practice worldwide. However, 3-D EM numerical simulations (modeling)—a key part of any 3-D EM data analysis—with realistic levels of complexity, accuracy and spatial detail still remains challenging from the computational point of view. We present a novel, efficient 3-D numerical solver based on a volume integral equation (IE) method. The efficiency is achieved by using a high-order polynomial (HOP) basis instead of the zero-order (piecewise constant) basis that is invoked in all routinely used IE-based solvers. We demonstrate that usage of the HOP basis allows us to decrease substantially the number of unknowns (preserving the same accuracy), with corresponding speed increase and memory saving. Electromagnetic theory, Geomagnetic induction, Magnetotellurics, Numerical modelling 1 INTRODUCTION In geophysics, electromagnetic (EM) methods aim to resolve generally 3-D distributions of subsurface electrical conductivity. Since conductivity depends on rock type and composition, temperature and fluid/melt content, these methods are widely used in academia and industry. To interpret EM data, the EM field which is governed by Maxwell’s equations needs to be computed for a given model of 3-D conductivity distribution. Here, we restrict our task by considering only the frequency-domain formulation. Usually, a large number of such simulations are required and complex, detailed and large-scale 3-D models are invoked. Thus, EM (forward) solvers that are able to deliver fast and accurate computation of the EM field are essential. The EM forward solvers can be distinguished by the numerical techniques utilized. Three basic techniques are used in geoelectromagnetism, namely finite-difference (FD), finite-element (FE) and volume integral equation (IE) methods. For decades FD-based solvers (Mackie et al.1994; Haber & Ascher 2001; Newman & Alumbaugh 2002; Egbert & Kelbert 2012, among others) dominated in EM due to the rather straightforward implementation of the FD concept. However, in recent years, FE-based (Schwarzbach et al.2011; Farquharson & Miensopust 2011; Puzyrev et al.2013; Ren et al.2013; Um et al.2013; Grayver & Burg 2014, among others) and IE-based (Avdeev et al.2002; Hursan & Zhdanov 2002; Singer 2008; Koyama et al.2008; Kamm & Pedersen 2014; Kruglyakov et al.2016, among others) solvers have increased in popularity due to methodological developments of the latter two methods. One of the main differences between the FE and IE methods lies in the structure and size of the resulting system matrices. In the IE method, one works with compact (but dense) matrices. The reason for compactness is that boundary conditions are exactly accounted for Green’s functions, and thus the modeling region is confined only to 3-D conductivity structures (anomalies) under investigations. By contrast, in the FE method, one has to discretize a much larger volume laterally and vertically in order to enable the decay (or stabilization) of the EM field at the boundaries of the modeling domain. However, this advantage is counterbalanced by the fact that FE matrices are sparse. Another distinction between the methods is that the condition number (which controls the stability of the solution) of system matrices in FE depends on discretization and frequency, whereas in IE—it does not, provided a particular class of IE with a contracting kernel is exploited (Singer 1995; Pankratov et al.1995; Zhdanov 2002). On the other hand, the FE method more easily treats the models with topography or/and bathymetry. It is pertinent to note here that all previous EM IE solvers (except the solver discussed in Farquharson & Oldenburg 2002), used the zero-order polynomial (piecewise constant; PWC) basis to describe the EM field behaviour within a modeling domain. Despite this, until recently, they demonstrated comparable efficiency with FE solvers based on a trilinear basis, in terms of accuracy and required computational resources. Recently, Grayver & Kolev (2015) showed that the usage of the second-order polynomial basis in FE solvers allows a decrease in the number of unknowns by several orders of magnitude and thus greatly speeds up EM field simulations. Inspired by this result, we developed an IE solver which exploits a high-order polynomial (HOP) basis, and demonstrate that it indeed outperforms the IE solver based on a PWC basis. 2 INTEGRAL EQUATION APPROACH Let the complex-valued function σ(M), |$\mathop {\mathbf {Re}}{\sigma (M)} \ge 0$| be a 3-D conductivity distribution in space. We search for the electric, E, and magnetic, H, fields induced by an electric current, Jext, in the model with conductivity σ(M). These fields obey the system of Maxwell’s equations \begin{eqnarray} \left\lbrace \begin{array}{ll} \mathop {\mathrm{curl}}\nolimits {\mathbf {H}}&= {\sigma }{\mathbf {E}}+{\mathbf {J}}_{\text{ext}},\\ \mathop {\mathrm{curl}}\nolimits {\mathbf {E}}&= i\omega \mu _0 {\mathbf {H}}.\\ \end{array} \right. \end{eqnarray} (1) Here, |$i=\sqrt{-1}$|, ω is an angular frequency, and μ0 is the magnetic permeability of free space. Time dependence of all fields is accounted for by e−iωt. Let |$\mathrm{\Omega }\subset \mathbb {R}^3$| be some bounded domain, σ(M) = σb(z) for |$M(x,y,z) \not\in \mathrm{\Omega }$|, and σ(M) = σa(M) for M ∈ Ω. Then, for any |$M \in \mathbb {R}^3$|, the fields E(M) and H(M) are expressed in terms of the integrals \begin{eqnarray} {\mathbf {E}}(M) ={\mathbf {E}}^b(M)+\int \limits _{\mathrm{\Omega }}\widehat{G}_E(M,M_0)\Delta _\sigma (M_0){\mathbf {E}}(M_0)d\mathrm{\Omega }_{M_0}, \end{eqnarray} (2) \begin{eqnarray} {\mathbf {H}}(M) ={\mathbf {H}}^b(M)+\int \limits _{\mathrm{\Omega }}\widehat{G}_H(M,M_0)\Delta _\sigma (M_0){\mathbf {E}}(M_0)d\mathrm{\Omega }_{M_0}. \end{eqnarray} (3) Here, Δσ = σa − σb and |$\widehat{G}_E$|, |$\widehat{G}_H$| are electric and magnetic Green’s tensors, respectively (cf. Pankratov et al.1995). The fields Eb and Hb in eqs (2) and (3) are called the background (normal) electric and magnetic fields and they obey the following system of Maxwell’s equations \begin{eqnarray} \left\lbrace \begin{array}{ll} &\mathop {\mathrm{curl}}\nolimits {\mathbf {H}}^b ={\sigma }_b(z) {\mathbf {E}}^b+{\mathbf {J}}_{\text{ext}},\\ &\mathop {\mathrm{curl}}\nolimits {\mathbf {E}}^b=i\omega \mu _0 {\mathbf {H}}^b . \\ \end{array} \right. \end{eqnarray} (4) Note, that the dependence on frequency of all the above quantities is omitted but implied. Note also, that |$\widehat{G}_E$| and |$\widehat{G}_H$| are dependent on the background conductivity σb. 2.1 Outline of Galerkin method Let us assume, that M ∈ Ω and rewrite eq. (2) in the following operator form \begin{eqnarray} (\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma ){\mathbf {E}}={\mathbf {E}}^b, \end{eqnarray} (5) where |$\mathop {\mathrm{{G}_{E}}}\nolimits$| is the integral operator from eq. (2) and |$\mathop {\mathrm{{I}_{}}}\nolimits$| is the identity operator. To solve eq. (5), we use the Galerkin method (Delves & Mohamed 1985; Farquharson & Oldenburg 2002; Zhdanov 2002) which is sketched below. Let |$\mathrm{\mathcal {L}_2}[\mathrm{\Omega }]$| be a vector Hilbert functional space with the following dot product and norm \begin{eqnarray} \left({\mathbf {V}},{\mathbf {U}}\right)&=& \int \nolimits _{\mathrm{\Omega }}\left(V_x(M)\overline{U}_x(M)+V_y(M)\overline{U}_y(M)\right.\nonumber\\ &&\left.+ \,V_z(M)\overline{U}_z(M)\right)d\mathrm{\Omega }_M,\nonumber\\ &&\left\Vert {\mathbf {V}}\right\Vert =\sqrt{\left({\mathbf {V}},{\mathbf {V}}\right)}. \end{eqnarray} (6) Here, the overline above U components denotes their complex conjugates. Let us also assume that vector functions Ψn, n = 1, …, N form the orthonormal basis in |$\mathrm{\mathcal {L}_2}[\mathrm{\Omega }]$|, denote |$\mathbf {\mathcal {Q}}^N$| as the linear span of this basis, and define projection operator |$\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits$| from |$\mathrm{\mathcal {L}_2}[\mathrm{\Omega }]$| to |$\mathbf {\mathcal {Q}}^N$| as follows \begin{eqnarray} {\mathbf {V}}^N&=&\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {V}}=\sum \limits _{n=1}^{N}a_n{\mathbf {\Psi }}_n,\nonumber\\ a_n &=&({\mathbf {V}},{\mathbf {\Psi }}_n),\nonumber\\ {\mathbf {V}} &\in& \mathrm{\mathcal {L}_2}[\mathrm{\Omega }] \quad {\mathbf {V}}^N \in \mathbf {\mathcal {Q}}^N. \end{eqnarray} (7) By construction |$\left\Vert \mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \right\Vert =1$| and VN is the best approximation of |${\mathbf {V}}\in \mathrm{\mathcal {L}_2}[\mathrm{\Omega }]$| in terms of norm (6). The Galerkin method as applied to eq. (5) relies on an approximation of E inside Ω by a function |${\mathbf {U}}\in \mathbf {\mathcal {Q}}^N$|, which satisfies the following N equations \begin{eqnarray} \left(\left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma \right){\mathbf {U}},{\mathbf {\Psi }}_n\right)=({\mathbf {E}}^b,{\mathbf {\Psi }}_n) \quad n=1,\ldots ,N. \end{eqnarray} (8) Substituting the expansion |${\mathbf {U}}=\sum \limits _{n=1}^{N}u_n{\mathbf {\Psi }}_n$| into eq. (8), we obtain the system of linear equations for coefficients un \begin{eqnarray} u_n-\sum \limits _{m=1}^{N}u_m(\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma {\mathbf {\Psi }}_m,{\mathbf {\Psi }}_n)=({\mathbf {E}}^b,{\mathbf {\Psi }}_n). \end{eqnarray} (9) By using the definition of projection operator |$\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits$|, we can deduce that function U obeys the following operator equation in space |$\mathbf {\mathcal {Q}}^N$| \begin{eqnarray} {\mathbf {U}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma {\mathbf {U}}=\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}^b. \end{eqnarray} (10) Equation (10) has a unique solution in |$\mathbf {\mathcal {Q}}^N$|, and thus the system (9) also has a unique solution (see more details in Appendix A). The final stage of the method is to substitute U instead of E into the integrands in eqs (2) and (3) in order to obtain |$\tilde{{\mathbf {E}}}$| and |$\tilde{{\mathbf {H}}}$| at any point |$M\in \mathbb {R}^3$| \begin{eqnarray} \tilde{{\mathbf {E}}}(M) ={\mathbf {E}}^b(M)+\int \nolimits _{\mathrm{\Omega }}\widehat{G}_E(M,M_0)\Delta _\sigma (M_0){\mathbf {U}}(M_0)\,d\mathrm{\Omega }_{M_0}, \end{eqnarray} (11) \begin{eqnarray} \tilde{{\mathbf {H}}}(M) ={\mathbf {H}}^b(M)+\int \nolimits _{\mathrm{\Omega }}\widehat{G}_H(M,M_0)\Delta _\sigma (M_0){\mathbf {U}}(M_0)\,d\mathrm{\Omega }_{M_0}. \end{eqnarray} (12) The approximations (11) and (12) guarantee |$\mathop {\mathrm{div}}\nolimits \tilde{{\mathbf {H}}}=0$| everywhere by design, irrespective of the choice of basis functions Ψn. Moreover, eqs (11) and (12) ensure that the components of |$\tilde{{\mathbf {E}}}$| and |$\tilde{{\mathbf {H}}}$| tangential to any plane are continuous across this plane. 2.2 Construction of basis functions It can be shown (see Appendix B) that for |${\mathbf {E}},\tilde{{\mathbf {E}}}$|, H and |$\tilde{{\mathbf {H}}}$| the following inequalities hold \begin{eqnarray} \left|\tilde{{\mathbf {E}}}(M)-{\mathbf {E}}(M)\right| \le C_E^{\text{out}}(M){\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert },\nonumber\\ \left|\tilde{{\mathbf {H}}}(M)-{\mathbf {H}}(M)\right| \le C_H^{\text{out}}(M){\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert }, \end{eqnarray} (13) for |$M \not\in \mathrm{\Omega }$| and \begin{eqnarray} {\left\Vert {\tilde{{\mathbf {E}}}-{\mathbf {E}}}\right\Vert } \le C_E^{\text{in}}{\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert },\nonumber\\ {\left\Vert {\tilde{{\mathbf {H}}}-{\mathbf {H}}}\right\Vert } \le C_H^{\text{in}}{\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert }, \end{eqnarray} (14) for M ∈ Ω. Here, |$C_{E,H}^{\text{out}}(M)$| are bounded functions and |$C_E^{\text{in}},C_H^{\text{in}}$| are constants, which are all independent of the basis functions Ψn. The inequalities (13) and (14) demonstrate rather the obvious but important fact that the accuracy of the numerical solution depends on |$\left\Vert {\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}\right\Vert$|, that is, how well the true electric field E is approximated by its projection to |$\mathbf {\mathcal {Q}}^N$|. An optimal choice of space |$\mathbf {\mathcal {Q}}^N$| is a non-trivial problem, but the insight gained in the FE modeling community (Schwarzbach et al.2011; Grayver & Kolev 2015) promotes the idea to use piecewise polynomial approximation of the electric field. The estimates above are based on properties of the solution, without considering efforts to calculate coefficients for system (9) quickly and accurately. Such a calculation is the main challenge in any IE approach, even for the PWC basis, where these coefficients are double volumetric integrals of Green’s tensors components (Kruglyakov & Bloshanskaya 2017). Note that the first attempt in the EM geophysical community to use the basis different from PWC was undertaken by Farquharson & Oldenburg (2002) who used trilinear approximation. However, the challenge to construct matrix coefficients for this basis compelled authors to use a Green’s tensor for the simplest background conductivity model, namely homogeneous space. In this paper, the calculation of the coefficients for the HOP basis is based on a generalization of the method proposed in Kruglyakov & Bloshanskaya (2017). Taking into account that the partial sum of Fourier–Legendre series provides the best polynomial approximation in |$\mathrm{\mathcal {L}_2}[-1,1]$|, we made an influential decision to use Legendre-type polynomials. Note, that in general, one can use basis functions from a wide class of Jacobi polynomials, but this implies the usage of weighting functions in dot products |$(\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma {\mathbf {\Psi }}_m,{\mathbf {\Psi }}_n)$|, which enormously complicates calculation of the matrix coefficients. The modeling domain Ω is split into Nc = NxNyNz non-overlapping rectangular cells |$\bigcup _{n=1}^{N_c} \mathrm{\Omega }_n=\mathrm{\Omega }$|; Nx,y,z is the number of cells in the x-, y- and z-directions, respectively, and conductivities σa,b are constant inside each cell. For every cell Ωn the finite scalar basis functions |$\Psi _n^{n_x,n_y,n_z}(x,y,z)$| are introduced as follows \begin{eqnarray} \Psi ^{n_x,n_y,n_z}_n(x,y,z)&=&\frac{2\sqrt{2}}{\sqrt{h_xh_yh_z}} L_{n_x}\left(2\frac{x-x^n}{h_x^n}-1\right)\nonumber\\ &&\times\, L_{n_y}\left(2\frac{y-y^n}{h_y^n}-1\right)\times L_{n_z}\left(2\frac{z-z^n}{h_z^n}-1\right)\nonumber\\ \end{eqnarray} (15) where, Lm is the normalized Legendre polynomial of the mth order, |$h_x^n$|, |$h_y^n$| and |$h_z^n$| are the sizes of Ωn in x, y and z dimensions, and xn, yn and zn are the coordinates of a corner of Ωn, |$n_{x,y,z}=0,\ldots ,N^P_{x,y,z}$|, and |$N^P_{x,y,z}$| is the maximum polynomial order along the x-, y- and z-directions. Note, that |$\Psi _n^{n_x,n_y,n_z}=0$| outside of Ωn, n = 1, …, Nc. The vector basis functions Ψn, n = 1, …, N, |$N=3N_c\left(N^P_x+1\right)\left(N^P_y+1\right)\left(N^P_z+1\right)$| are constructed as Ψ = (Ψx, Ψy, Ψz), where Ψx, y, z are those of |$\Psi ^{n_x,n_y,n_z}_{n}$|. The orthonormality of Ψn follows from the orthonormality of (normalized) Legendre polynomials. In the terminology of FE, we use the discontinuous Galerkin method, with basis functions that do not guarantee continuity of U at Ωn boundaries. However, the continuity of tangential components of |$\tilde{{\mathbf {E}}}$| and |$\tilde{{\mathbf {H}}}$| across the cells’ boundaries is secured by eqs (11) and (12), as well as zero divergence of |$\tilde{{\mathbf {H}}}$|. 3 CONVERGENCE RATE OF THE METHOD The inequalities (13) and (14) are formulated in terms of |${\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert }$|. We can make a step further and estimate |${\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert }$| in terms of the maximum cell diameter d and the minimum polynomial order |$N^P=\min \lbrace N_x^P,N^P_y,N^P_z\rbrace$| \begin{eqnarray} {\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert }_{\mathrm{\mathcal {L}_2}[\mathrm{\Omega }]} \le A d^{N^P+1}, \end{eqnarray} (16) where A is a constant independent of d (but dependent on NP and E). This result is a corollary from the Bramble–Hilbert lemma (Brenner & Scott 2008), if one recalls that |$\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}$| is the best polynomial approximation in |$\mathrm{\mathcal {L}_2}[\mathrm{\Omega }_n]$| by design. Substituting inequality (16) into the right-hand side of inequalities (13) and (14), we obtain that for the fixed polynomial order NP the method converges with NP + 1 order. 4 REMARKS ON IMPLEMENTATION The main challenges of any IE-based solver are computation of matrix coefficients [see eq. (9)] and (dense) matrix handling. These challenges have been addressed by using the approaches introduced by Kruglyakov & Bloshanskaya (2017). More specifically, the integrals in the vertical direction are computed analytically by generalization of the approach of Kruglyakov & Bloshanskaya (2017, appendix B), whereas for the horizontal integration special digital filters are constructed (appendix C of the same paper). However, the analytical transformations involved at these stages are extremely complex and we used the computer algebra system Maxima (2017) to obtain the correct expressions and to generate proper Fortran code. A standard IE numerical scheme is invoked, based on a uniform discretization in lateral directions. It yields a block-Toeplitz system matrix, thus decreasing memory requirements and accelerating matrix-vector multiplication by using a 2-D fast Fourier transform (cf. Avdeev et al.1997). In addition, system matrix symmetries and antisymmetries were utilized to further reduce computational and memory loads as in Kruglyakov & Bloshanskaya (2017). Finally, the resulting system of equations is solved by Krylov subspace iterations using our own implementation of the FGMRES method (Saad 1993). To obtain faster convergence to the solution, we use IE with a contracting kernel approach (Singer 1995; Pankratov et al.1995; Zhdanov 2002), hereinafter denoted as contracting IE (CIE). The advantage of using CIE is that the condition number of the CIE system matrix depends only on the lateral conductivity contrast (cf. Pankratov & Kuvshinov 2016), which permits use of the HOP method without additional preconditioning of the system matrix. The only reason we use conventional IE in the equations above is to simplify the explanation. 5 NUMERICAL EXPERIMENTS To explore the performance of the HOP IE-based solver against the one which uses PWC basis functions, we performed simulations of magnetotelluric (MT) responses using three different models of 3-D conductivity distribution. 5.1 Fault model We start with presenting results for a case in which polynomials are implemented along the vertical direction only. The model taken from Bakker et al. (2015) mimics a conducting sedimentary basin surrounded by resistive mountains with a conductive fault beneath the left flank of the basin. Plane and side views of the model are shown in the respective upper and lower plots in Fig. 1. We expect galvanic coupling between the high-conductive sedimentary basin and the deep conducting basement by the downward leakage of electric currents within the fault, thus assuming that accurate approximation of the electric field in the vertical direction is essential. Figure 1. Open in new tabDownload slide Fault model. Top (a)—a plan view of the model and bottom (b)—side view. The dashed line in the upper plot depicts the profile along which MT apparent resistivities are shown in Fig. 2. Figure 1. Open in new tabDownload slide Fault model. Top (a)—a plan view of the model and bottom (b)—side view. The dashed line in the upper plot depicts the profile along which MT apparent resistivities are shown in Fig. 2. The need for accurate approximation is clearly seen in Fig. 2, where the modeling by a PWC IE solver with vertical discretization Nz = 32 produces an artefact at the right flank of the fault. Increasing Nz in the course of PWC modeling improves the results, as expected. The artefact disappears when rather excessive discretization with Nz = 128 is invoked. Remarkably, the usage of a third-order polynomial basis in the vertical direction gives accurate results already with Nz = 2. Note that lateral discretization hx = hy = 200 m was used for all runs, for both PWC and HOP modeling. Details of the successful runs are summarized in Table 1. It is evident from the table that in order to reach the same accuracy one needs 16 times fewer unknowns if the HOP solver is invoked. In terms of computational loads the gains are 25 and 8 for the memory and time, respectively. Figure 2. Open in new tabDownload slide Apparent resistivity ρxy (at period 1000 s) along the profile shown in Fig. 1(a). Solid (blue, red and yellow) curves are the results obtained with the use of the PWC IE solver. The colours distinguish between the results of modelings with different numbers of cells in the vertical direction, Nz. The dashed curve stands for the results obtained with the use of the HOP IE solver. Figure 2. Open in new tabDownload slide Apparent resistivity ρxy (at period 1000 s) along the profile shown in Fig. 1(a). Solid (blue, red and yellow) curves are the results obtained with the use of the PWC IE solver. The colours distinguish between the results of modelings with different numbers of cells in the vertical direction, Nz. The dashed curve stands for the results obtained with the use of the HOP IE solver. Table 1. Details of successful runs with HOP and PWC IE solvers for the model depicted in Fig. 1. |$\tilde{N}_z$| stands for the number of unknowns in vertical direction. In PWC case |$\tilde{N}_z = N_z$|, in HOP case |$\tilde{N}_z = 4N_z$|, where 4 arises since the third-order polynomial basis is used. Time is that required for two polarizations modeling. Run with |$\tilde{N}_z$| Number of unknowns RAM (GB) Time (CPU × hours) PWC 128 96 × 106 427 17 HOP 8 6 × 106 17 2 Run with |$\tilde{N}_z$| Number of unknowns RAM (GB) Time (CPU × hours) PWC 128 96 × 106 427 17 HOP 8 6 × 106 17 2 Open in new tab Table 1. Details of successful runs with HOP and PWC IE solvers for the model depicted in Fig. 1. |$\tilde{N}_z$| stands for the number of unknowns in vertical direction. In PWC case |$\tilde{N}_z = N_z$|, in HOP case |$\tilde{N}_z = 4N_z$|, where 4 arises since the third-order polynomial basis is used. Time is that required for two polarizations modeling. Run with |$\tilde{N}_z$| Number of unknowns RAM (GB) Time (CPU × hours) PWC 128 96 × 106 427 17 HOP 8 6 × 106 17 2 Run with |$\tilde{N}_z$| Number of unknowns RAM (GB) Time (CPU × hours) PWC 128 96 × 106 427 17 HOP 8 6 × 106 17 2 Open in new tab 5.2 Two block (COMMEMI3D-2) model The developed HOP IE solver permits the use of polynomials up to fifth order in all directions. In this section, we explore the interplay between the order of polynomials in a lateral direction and lateral discretization of the model. Modelings were performed using the COMMEMI3D-2 model from Zhdanov et al. (1997) (see Fig. 3) in the period range between 1 and 100 s. In all HOP simulations, we use third-order polynomials in the vertical direction and Nz = 4 with vertical cell sizes of 500, 500, 4000 and 5000 m. In the course of HOP simulations, the following combinations of lateral discretization of the model, NH, and order of polynomials in the lateral direction, |$N^P_H$|, were tried: |$N_H = 4, N_H^P=5$|; |$N_H = 8, N_H^P=5$| and |$N_H = 16, N_H^P=3$|. For all PWC simulations, NH = 128 and Nz = 50 with vertical cell sizes growing in geometric sequence from 10 to 845 m. Numerical experiments (not shown here) demonstrate that this combination of NH and Nz provides accurate PWC results for all considered periods. Figure 3. Open in new tabDownload slide COMMEMI3D-2 model. Top (a)—a plan view of the model and bottom (b)—side view. Dashed line in the top plot depicts the profile along which MT responses are shown in Figs 4 and 5. Figure 3. Open in new tabDownload slide COMMEMI3D-2 model. Top (a)—a plan view of the model and bottom (b)—side view. Dashed line in the top plot depicts the profile along which MT responses are shown in Figs 4 and 5. The numerical simulations also show that the efficiency of HOP may depend on the period. At a period of 100 s (Fig. 4), the electric field inside a cell is rather smooth and the HOP approach yields accurate results using a very coarse lateral grid (two cells per block) provided fifth-order polynomials are used. At the same time, the PWC solver requires 128 cells in a lateral direction to obtain results of comparable accuracy. Cumulatively, for a period of 100 s, by using the HOP solver one needs 100 times fewer unknowns to obtain accurate results, and the corresponding gains in terms of computational loads are 6 and 7 for the memory and time, respectively. Figure 4. Open in new tabDownload slide Apparent resistivities and phases along the profile shown by the dashed line in Fig. 3(a). Period is 100 s. Figure 4. Open in new tabDownload slide Apparent resistivities and phases along the profile shown by the dashed line in Fig. 3(a). Period is 100 s. At a period of 1 s, accurate results by the HOP are obtained for the combination |$N_H = 16, N_H^P=3$| meaning that a denser lateral grid is used (Fig. 5). At this period, the electric field varies sharply inside the conductive block (Fig. 6) and thus the lateral cell’s size has to be decreased (which is what we did) or the polynomial order has to be significantly increased. Figure 5. Open in new tabDownload slide The same legend as in Fig. 4. Period is 1 s. Figure 5. Open in new tabDownload slide The same legend as in Fig. 4. Period is 1 s. Figure 6. Open in new tabDownload slide Magnitude of Ex for x-polarization excitation at z = 0 along the discussed profile over the left half of the conductive block in the COMMEMI3D-2 model at a period of 1 s. Figure 6. Open in new tabDownload slide Magnitude of Ex for x-polarization excitation at z = 0 along the discussed profile over the left half of the conductive block in the COMMEMI3D-2 model at a period of 1 s. Despite a denser lateral grid for this period (NH = 16), the HOP solver still outperforms PWC solver, but the gain is smaller (3 and 1.8 for the memory and time) than for 100 s modeling. However, we note here that the HOP performance can be improved by using non-regular lateral grids, which is the subject of further research. Details of successful HOP and PWC modeling runs are summarized in Table 2. Table 2. Details of successful modeling runs with HOP and PWC IE solvers for the COMMEMI3D-2 model. NH stands for a number of cells in every lateral direction. |$N_H^P$| is the polynomial order used in the lateral direction. Time is that required for modeling two polarizations. Run with NH |$N_H^P$| Number of unknowns RAM (GB) Time (CPU × s) HOP (100 s) 4 5 27 648 1.5 195 HOP (1 s) 16 3 196 608 3.5 714 PWC 128 2457 600 10 1260 Run with NH |$N_H^P$| Number of unknowns RAM (GB) Time (CPU × s) HOP (100 s) 4 5 27 648 1.5 195 HOP (1 s) 16 3 196 608 3.5 714 PWC 128 2457 600 10 1260 Open in new tab Table 2. Details of successful modeling runs with HOP and PWC IE solvers for the COMMEMI3D-2 model. NH stands for a number of cells in every lateral direction. |$N_H^P$| is the polynomial order used in the lateral direction. Time is that required for modeling two polarizations. Run with NH |$N_H^P$| Number of unknowns RAM (GB) Time (CPU × s) HOP (100 s) 4 5 27 648 1.5 195 HOP (1 s) 16 3 196 608 3.5 714 PWC 128 2457 600 10 1260 Run with NH |$N_H^P$| Number of unknowns RAM (GB) Time (CPU × s) HOP (100 s) 4 5 27 648 1.5 195 HOP (1 s) 16 3 196 608 3.5 714 PWC 128 2457 600 10 1260 Open in new tab 5.3 Large conductivity contrast (COMMEMI3D-3) model An accurate computation of electric and magnetic fields in media with a large conductivity contrast is one of the most challenging problems of EM modeling. This is due to strong interconnection between the lateral conductivity contrast and the matrix condition number (Pankratov & Kuvshinov 2016). The large conductivity contrast COMMEMI3D-3 model (Hursan & Zhdanov 2002; Varentsov et al.2000) is considered as the third test model for the presented solver. This model schematically describes the conductivity distribution typical for ore exploration with MT. The COMMEMI3D-3 model consists of seven rectangular blocks embedded in a layered medium (Fig. 7). The background section of the model consists of two layers with conductivities of 10−3, 10−4 S m−1, and a lower half-space with conductivity of 0.1 S m−1. The thicknesses of the first and the second layers are 1 and 6.5 km, respectively. The maximum lateral conductivity contrast is 104 in the first layer and 3.3 × 104 in the second layer. Figure 7. Open in new tabDownload slide COMMEMI3D-3 model. Left-hand plots show side views of the model and right-hand plots show the plane view. Figure 7. Open in new tabDownload slide COMMEMI3D-3 model. Left-hand plots show side views of the model and right-hand plots show the plane view. The MT responses were calculated along two representative profiles (Fig. 7) in a period range between 0.01 and 1000 s using different discretizations and polynomial orders (cf. Table 3). The IE modeled responses (Figs 8–13) were compared with those from the reference second-order FE solver by Grayver & Kolev (2015). Figure 8. Open in new tabDownload slide Apparent resistivities ρxy, ρyx (above) and impedance phases φxy, φyx (below) along the profile y = 3830. Period is 1 s. Figure 8. Open in new tabDownload slide Apparent resistivities ρxy, ρyx (above) and impedance phases φxy, φyx (below) along the profile y = 3830. Period is 1 s. Figure 9. Open in new tabDownload slide Real (left) and imaginary (right) parts of tippers component Wxz, Wyz along the profile y = 3830. Period is 1 s. Figure 9. Open in new tabDownload slide Real (left) and imaginary (right) parts of tippers component Wxz, Wyz along the profile y = 3830. Period is 1 s. Figure 10. Open in new tabDownload slide Apparent resistivities ρxy, ρyx (above) and impedance phases φxy, φyx (below) along the profile x = 1900. Period is 1 s. Figure 10. Open in new tabDownload slide Apparent resistivities ρxy, ρyx (above) and impedance phases φxy, φyx (below) along the profile x = 1900. Period is 1 s. Figure 11. Open in new tabDownload slide Real (left) and imaginary (right) parts of tippers component Wxz, Wyz along the profile x = 1900. Period is 1 s. Figure 11. Open in new tabDownload slide Real (left) and imaginary (right) parts of tippers component Wxz, Wyz along the profile x = 1900. Period is 1 s. Figure 12. Open in new tabDownload slide Apparent resistivities ρxy, ρyx (above) and impedance phases φxy, φyx (below) with respect to period at point (3975, 3830). Figure 12. Open in new tabDownload slide Apparent resistivities ρxy, ρyx (above) and impedance phases φxy, φyx (below) with respect to period at point (3975, 3830). Figure 13. Open in new tabDownload slide Real (left) and imaginary (right) parts of tippers component Wxz, Wyz (below) with respect to period at point (3975, 3830). Figure 13. Open in new tabDownload slide Real (left) and imaginary (right) parts of tippers component Wxz, Wyz (below) with respect to period at point (3975, 3830). Table 3. Details of modeling runs for COMMEMI3D-3 model with HOP and PWC IE solvers. |$N_H^P$| and |$N_z^P$| stand for polynomial orders used in lateral and vertical directions, respectively. Legend Run with Discretization |$N_H^P$| |$N_z^P$| Number of unknowns PWC coarse PWC hx = hy = 100 m, hz from 25 to 250 m 177 408 PWC fine PWC hx = hy = hz = 12.5 m 111648 768 PWC best PWC hx = hy = hz = 6.25 m 893190 144 HOP |$N_H^P=1$| HOP hx = hy = 200 m, hz from 102 to 103 m 1 3 177 408 HOP |$N_H^P=2$| HOP hx = hy = 200 m, hz from 102 to 103 m 2 3 399 168 Legend Run with Discretization |$N_H^P$| |$N_z^P$| Number of unknowns PWC coarse PWC hx = hy = 100 m, hz from 25 to 250 m 177 408 PWC fine PWC hx = hy = hz = 12.5 m 111648 768 PWC best PWC hx = hy = hz = 6.25 m 893190 144 HOP |$N_H^P=1$| HOP hx = hy = 200 m, hz from 102 to 103 m 1 3 177 408 HOP |$N_H^P=2$| HOP hx = hy = 200 m, hz from 102 to 103 m 2 3 399 168 Open in new tab Table 3. Details of modeling runs for COMMEMI3D-3 model with HOP and PWC IE solvers. |$N_H^P$| and |$N_z^P$| stand for polynomial orders used in lateral and vertical directions, respectively. Legend Run with Discretization |$N_H^P$| |$N_z^P$| Number of unknowns PWC coarse PWC hx = hy = 100 m, hz from 25 to 250 m 177 408 PWC fine PWC hx = hy = hz = 12.5 m 111648 768 PWC best PWC hx = hy = hz = 6.25 m 893190 144 HOP |$N_H^P=1$| HOP hx = hy = 200 m, hz from 102 to 103 m 1 3 177 408 HOP |$N_H^P=2$| HOP hx = hy = 200 m, hz from 102 to 103 m 2 3 399 168 Legend Run with Discretization |$N_H^P$| |$N_z^P$| Number of unknowns PWC coarse PWC hx = hy = 100 m, hz from 25 to 250 m 177 408 PWC fine PWC hx = hy = hz = 12.5 m 111648 768 PWC best PWC hx = hy = hz = 6.25 m 893190 144 HOP |$N_H^P=1$| HOP hx = hy = 200 m, hz from 102 to 103 m 1 3 177 408 HOP |$N_H^P=2$| HOP hx = hy = 200 m, hz from 102 to 103 m 2 3 399 168 Open in new tab In all figures, one observes excellent agreement between the HOP IE and the FE solvers. Note that this agreement is achieved with the HOP IE solver by using a very moderate number (around 2 × 105) of unknowns, with bilinear basis in lateral directions and cubic basis in the vertical direction. Note that the increasing polynomial order in lateral directions practically does not change the results. Interestingly, when the HOP solver is used, a rather coarse lateral grid (with hx = hy = 200 m) is sufficient for accurate modeling of the responses over the shallow structure. It is quite unexpected, because the electric field varies non-negligibly inside shallow blocks and their sizes in the y-direction (from 400 to 600 m) are comparable with the cell sizes (Figs 10–11). At the same time, using the PWC IE solver requires around 9 × 108 unknowns to obtain results of comparable accuracy, which leads to a corresponding difference in computational loads (cf. Table 4). A coarser discretization does not produce ‘true’ results in apparent resistivity over the conductive block at periods from 0.1 to 10 s (upper parts of Figs 8 and 12). Nor does it do so for the phases at any periods (lower parts of Figs 8 and 12). Moreover, even with the ‘best’ discretization (cubic cells of 6.25 m size), the PWC impedance phases do not perfectly agree with those from the FE and HOP solvers. As a final remark, if PWC and HOP discretizations are comparable (by giving exactly the same number of unknowns), the PWC results differ considerably from the ‘true’ results, especially above the high-conductive block, and most prominently in phases of impedance and tippers (Figs 8–13). Table 4. Details of successful modeling runs for COMMEMI3D-3 model with HOP and PWC IE solvers. Time is for modeling two polarizations. Legend Number of unknowns RAM (GB) Time (CPU × hours) PWC best 893190 144 10 000 4500 HOP |$N_H^P=1$| 177 408 2.81 0.8 Legend Number of unknowns RAM (GB) Time (CPU × hours) PWC best 893190 144 10 000 4500 HOP |$N_H^P=1$| 177 408 2.81 0.8 Open in new tab Table 4. Details of successful modeling runs for COMMEMI3D-3 model with HOP and PWC IE solvers. Time is for modeling two polarizations. Legend Number of unknowns RAM (GB) Time (CPU × hours) PWC best 893190 144 10 000 4500 HOP |$N_H^P=1$| 177 408 2.81 0.8 Legend Number of unknowns RAM (GB) Time (CPU × hours) PWC best 893190 144 10 000 4500 HOP |$N_H^P=1$| 177 408 2.81 0.8 Open in new tab 6 CONCLUDING REMARKS We present an IE solver which exploits high- (up to fifth) order polynomial basis functions. The solver outperforms IE solvers based on a conventional PWC basis in terms of both memory savings and acceleration of computations due to a decrease in the number of unknowns by one to several orders of magnitude. The gain depends on model and period. In the best scenario presented in the paper, the gain is four orders of magnitude. Such extreme reduction in computational loads opens an avenue for using non-uniform grids in lateral directions which were computationally forbidden in PWC IE solvers. It will allow, in particular, to use the IE approach for large-scale regional modeling where complex local 3-D effects from topography and bathymetry are important. We applied the solver to an MT problem. The extension of the code on the controlled-source case is rather straightforward and is the topic of ongoing realization. Finally, we note that the presented solver is designed to work in Cartesian geometry. The extension of the concept to a spherical geometry is also the topic of future research. Acknowledgements The authors thank Alexander Grayver for providing his FE results and William Lowrie for improving the English presentation of this paper. We also extend our gratitude to Michael Zhdanov, Takao Koyama and Joerg Renner for many constructive comments on the manuscript. This work is supported by the Swiss Government Excellence Scholarships for Foreign Scholars, Scholarship no 2016.0049 and the Swiss National Supercomputing Center (CSCS) grant (project ID s767). REFERENCES Avdeev D. , Kuvshinov A. , Pankratov O. , Newman G. , 1997 . High-performance three-dimensional electromagnetic modeling using modified Neumann series. Wide-band numerical solution and examples , J. Geomagn. Geoelectr. , 49 , 1519 – 1539 . https://doi.org/10.5636/jgg.49.1519 Google Scholar Crossref Search ADS WorldCat Avdeev D. , Kuvshinov A. , Pankratov O. , Newman G. , 2002 . Three-dimensional induction logging problems, Part I: an integral equation solution and model comparisons , 67 ( 2 ), 413 – 426 . https://doi.org/10.1190/1.1468601 WorldCat Bakker J. , Kuvshinov A. , Samrock F. , Geraskin A. , Pankratov O. , 2015 . Introducing inter-site phase tensors to suppress galvanic distortion in the telluric method , Earth Planets Space , 67 ( 1 ), 160 , doi:10.1186/s40623-015-0327-7 . https://doi.org/10.1186/s40623-015-0327-7 Google Scholar Crossref Search ADS WorldCat Brenner S. , Scott L. , 2008 . The Mathematical Theory of Finite Element Methods , Springer New York . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Delves L. , Mohamed J. , 1985 . Computational Methods for Integral Equations , Cambridge University Press . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Egbert G.D. , Kelbert A. , 2012 . Computational recipes for electromagnetic inverse problems , Geophys. J. Int. , 189 ( 1 ), 251 – 267 . https://doi.org/10.1111/j.1365-246X.2011.05347.x Google Scholar Crossref Search ADS WorldCat Farquharson C.G. , Miensopust M.P. , 2011 . Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction , J. appl. Geophys. , 75 , 699 – 710 . https://doi.org/10.1016/j.jappgeo.2011.09.025 Google Scholar Crossref Search ADS WorldCat Farquharson C.G. , Oldenburg D.W. , 2002 . An integral equation solution to the geophysical electromagnetic forward-modelling problem , in Three-Dimensional Electromagnetics: Methods in Geochemistry and , 35 , pp. 3 – 19 , eds Zhdanov M.S. , Wannamaker P.E. , Elsevier . WorldCat Gilbarg D. , Trudinger N. , 2001 . Elliptic Partial Differential Equations of Second Order , Springer , Berlin Heidelberg . Google Preview WorldCat COPAC Grayver A. , Burg M. , 2014 . Robust and scalable 3-D geo-electromagnetic modelling approach using the finite element method , Geophys. J. Int. , 198 , 110 – 125 . https://doi.org/10.1093/gji/ggu119 Google Scholar Crossref Search ADS WorldCat Grayver A. , Kolev T. , 2015 . Large-scale 3D geo-electromagnetic modeling using parallel adaptive high-order finite element method , 80 ( 6 ), 277 – 291 . https://doi.org/10.1190/geo2015-0013.1 WorldCat Haber E. , Ascher U.M. , 2001 . Fast finite volume simulation of 3D electromagnetic problems with highly discontinuous coefficients , SIAM J. Sci. Comput. , 22 ( 6 ), 1943 – 1961 . https://doi.org/10.1137/S1064827599360741 Google Scholar Crossref Search ADS WorldCat Hursan G. , Zhdanov M.S. , 2002 . Contraction integral equation method in three-dimensional electromagnetic modeling , Radio Sci. , 37 ( 6 ), 1–13, 1089. https://doi.org/10.1029/2001RS002513 WorldCat Kamm J. , Pedersen L.B. , 2014 . Inversion of airborne tensor VLF data using integral equations , Geophys. J. Int. , 198 ( 2 ), 775 – 794 . https://doi.org/10.1093/gji/ggu161 Google Scholar Crossref Search ADS WorldCat Koyama T. , Utada H. , Avdeev D. , 2008 . Fast and memory-saved 3-D forward modeling code for MT by using integral equation method , in Abstract book , p. 725 , ed. Chinese Geophysical Society, 19th Workshop on Electromagnetic Induction in the Earth , Beijing,China . Google Preview WorldCat COPAC Kruglyakov M. , Bloshanskaya L. , 2017 . High-performance parallel solver for integral equations of electromagnetics based on Galerkin method , Math. Geosci. , 49 ( 6 ), 751 – 776 . https://doi.org/10.1007/s11004-017-9677-y Google Scholar Crossref Search ADS WorldCat Kruglyakov M. , Geraskin A. , Kuvshinov A. , 2016 . Novel accurate and scalable 3-D MT forward solver based on a contracting integral equation method , Comput. Geosci. , 96 , 208 – 217 . https://doi.org/10.1016/j.cageo.2016.08.017 Google Scholar Crossref Search ADS WorldCat Mackie R. , Smith J. , Madden T. , 1994 . 3-Dimensional electromagnetic modeling using finite-difference equation—the magnetotelluric example , Radio Sci. , 29 ( 4 ), 923 – 935 . https://doi.org/10.1029/94RS00326 Google Scholar Crossref Search ADS WorldCat Maxima , 2017 . Maxima, a Computer Algebra System. Version 5.40 , http://maxima.sourceforge.net/, last access date 28 February 2018 . Google Preview WorldCat COPAC Newman G. , Alumbaugh D. , 2002 . Three-dimensional induction logging problems, Part 2: a finite-difference solution , 61 , 484 – 491 . https://doi.org/10.1190/1.1468608 WorldCat Pankratov O. , Kuvshinov A. , 2016 . Applied mathematics in EM studies with special emphasis on an uncertainty quantification and 3-D integral equation modelling , Surv. Geophys. , 37 ( 1 ), 109 – 147 . https://doi.org/10.1007/s10712-015-9340-4 Google Scholar Crossref Search ADS WorldCat Pankratov O. , Avdeyev D. , Kuvshinov A. , 1995 . Electromagnetic-field scattering in a heterogeneous Earth: a solution to the forward problem , Izv-Phys. Solid Earth , 31 ( 3 ), 201 – 209 . WorldCat Puzyrev V. , Koldan J. , de la Puente J. , Houzeaux G. , Vazquez M. , Cela J.M. , 2013 . A parallel finite-element method for three-dimensional controlled-source electromagnetic forward modelling , Geophys. J. Int. , 193 , 678 – 693 . https://doi.org/10.1093/gji/ggt027 Google Scholar Crossref Search ADS WorldCat 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 , 700 – 718 . https://doi.org/10.1093/gji/ggt154 Google Scholar Crossref Search ADS WorldCat Saad Y. , 1993 . A flexible inner-outer preconditioned GMRES algorithm , SIAM J. Sci. Comput. , 14 ( 2 ), 461 – 469 . https://doi.org/10.1137/0914028 Google Scholar Crossref Search ADS WorldCat 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 , 63 – 74 . https://doi.org/10.1111/j.1365-246X.2011.05127.x Google Scholar Crossref Search ADS WorldCat Singer B. , 1995 . Method for solution of Maxwell’s equations in non-uniform media , Geophys J Int , 120 , 590 – 598 . https://doi.org/10.1111/j.1365-246X.1995.tb01841.x Google Scholar Crossref Search ADS WorldCat Singer B. , 2008 . Electromagnetic integral equation approach based on contraction operator and solution optimization in Krylov subspace , Geophys. J. Int. , 175 , 857 – 884 . https://doi.org/10.1111/j.1365-246X.2008.03930.x Google Scholar Crossref Search ADS WorldCat Um E.S. , Commer M. , Newman G.A. , 2013 . Efficient pre-conditioned iterative solution strategies for the electromagnetic diffusion in the earth: finite-element frequency-domain approach , Geophys. J. Int. , 193 , 1460 – 1473 . https://doi.org/10.1093/gji/ggt071 Google Scholar Crossref Search ADS WorldCat Varentsov I.M. , Fomenko I.Y. , Golubev N.G. , Mehanee S. , Hursan G. , Zhdanov M.S. , 2000 . Comparative study of 3-D finite difference and integral equation methods , in Proceedings of 2000 Consortium for Electomagnetic Modeling and Inversion Annual Meeting , pp. 35 – 74 , University of Utah , Salt Lake City . Google Preview WorldCat COPAC Ward S.H. , Hohmann G.W. , 1988 . 4. Electromagnetic Theory for Geophysical Applications , chap. 4 , pp. 130 – 311 , SEG , USA . Google Preview WorldCat COPAC Zhdanov M. , Varentsov I. , Weaver J. , Golubev N. , Krylov V. , 1997 . Methods for modelling electromagnetic fields results from commemi-the international project on the comparison of modelling methods for electromagnetic induction , J. appl. Geophys. , 37 ( 3 ), 133 – 271 . https://doi.org/10.1016/S0926-9851(97)00013-X Google Scholar Crossref Search ADS WorldCat Zhdanov M.S. , 2002 . Chapter 9 - integral representations in electromagnetic forward modeling , in Geophysical Inverse Theory and Regularization Problems, Methods in Geochemistry and , 36 , pp. 231 – 286 , ed. Zhdanov M.S. , Elsevier . Google Preview WorldCat COPAC APPENDIX A: EXISTENCE AND UNIQUENESS OF THE IE SOLUTION BASED ON GALERKIN METHOD The simplest way to prove the existence and uniqueness of the solution of eq. (10) is to use the CIE approach. Let us make the natural assumption that |$\mathop {\mathbf {Re}}{{\sigma }_b}(z)>0$| for M(x, y, z) ∈ Ω and define operator |$\mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}$| as follows \begin{eqnarray} \mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}\mathbf { V}_{{}}=\sqrt{\mathop {\mathbf {Re}}{{\sigma }_b}}\mathop {\mathrm{{G}_{E}}}\nolimits \left[2\sqrt{\mathop {\mathbf {Re}}{{\sigma }_b}}\mathbf { V}_{{}}\right]+\mathbf { V}_{{}}. \end{eqnarray} (A1) Using eq. (A1), we can express |$\mathop {\mathrm{{G}_{E}}}\nolimits$| as \begin{eqnarray} \mathop {\mathrm{{G}_{E}}}\nolimits =\frac{1}{2\sqrt{\mathop {\mathbf {Re}}{{\sigma }_b}}}\left(\mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}-\mathop {\mathrm{{I}_{}}}\nolimits \right)\frac{1}{\sqrt{\mathop {\mathbf {Re}}{{\sigma }_b}}} \end{eqnarray} (A2) and rewrite eq. (10) as \begin{eqnarray} \left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}\frac{b}{a}\right)\tilde{{\mathbf {U}}}=\sqrt{\mathop {\mathbf {Re}}{{\sigma }_b}}{\mathbf {E}}^N, \end{eqnarray} (A3) where \begin{eqnarray} \tilde{{\mathbf {U}}}=a{\mathbf {U}},\quad a=\frac{{\sigma }_a+\overline{{\sigma }_b}}{2\sqrt{\mathop {\mathbf {Re}}{{\sigma }_b}}}\quad b=\frac{{\sigma }_a-{\sigma }_b}{2\sqrt{\mathop {\mathbf {Re}}{{\sigma }_b}}}. \end{eqnarray} (A4) Taking into account that |${\left\Vert {\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits }\right\Vert }=1,\left|\frac{b}{a}\right|<1$|, and |$\mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}$| is a contracting operator (Pankratov et al.1995; Singer 1995), we can construct |$\left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}\frac{b}{a}\right)^{-1}$| as the following Neumann’s series \begin{eqnarray} \left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}\frac{b}{a}\right)^{-1}=\sum \limits _{n=0}^{\infty }\left(\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}\frac{b}{a}\right)^n. \end{eqnarray} (A5) Using eq. (A5), we obtain \begin{eqnarray} {\left\Vert {\left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}\frac{b}{a}\right)^{-1}}\right\Vert } \le \sum \limits _{n=0}^{\infty }{\left\Vert {\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}\frac{b}{a}}\right\Vert }^n. \end{eqnarray} (A6) Hence, using definitions (A4), and taking into account that |$\left|\frac{{\sigma }_a-{\sigma }_b}{{\sigma }_a+\overline{{\sigma }_b}}\right|<1$| [cf. Pankratov & Kuvshinov (2016)] gives \begin{eqnarray} \sum \limits _{n=0}^{\infty }{\left\Vert {\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}\frac{b}{a}}\right\Vert }^n \le \sum \limits _{n=0}^{\infty }\left(\max \limits _{\mathrm{\Omega }}\left|\frac{{\sigma }_a-{\sigma }_b}{{\sigma }_a+\overline{{\sigma }_b}}\right|\right)^n=\frac{1}{1- \max \limits _{\mathrm{\Omega }}\left|\frac{{\sigma }_a-{\sigma }_b}{{\sigma }_a+\overline{{\sigma }_b}}\right|}.\nonumber\\ \end{eqnarray} (A7) Using eq. (A1) to express |$\mathop {\mathrm{{G}_{E}}}\nolimits$| in terms of |$\mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}$|, we have \begin{eqnarray} \left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma \right)^{-1}=2\frac{\sqrt{\mathop {\mathbf {Re}}{{\sigma }_b}}}{{\sigma }_a+\overline{{\sigma }_b}}\left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{\mathbf {G}_{E}^{m}}}}\nolimits {}\frac{b}{a}\right)^{-1}\sqrt{\mathop {\mathbf {Re}}{{\sigma }_b}}. \end{eqnarray} (A8) Thus, we conclude that operator |$\left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma \right)^{-1}$| is bounded \begin{eqnarray} {\left\Vert {\left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma \right)^{-1}}\right\Vert } \le S=2\max \limits _{\mathrm{\Omega }}\frac{\sqrt{\mathop {\mathbf {Re}}{{\sigma }_b}}}{\left|{\sigma }_a+\overline{{\sigma }_b}\right|}\frac{\max \limits _{\mathrm{\Omega }}\sqrt{\mathop {\mathbf {Re}}{{\sigma }_b}}}{1- \max \limits _{\mathrm{\Omega }}\left|\frac{{\sigma }_a-{\sigma }_b}{{\sigma }_a+\overline{{\sigma }_b}}\right|}. \nonumber\\ \end{eqnarray} (A9) and consequently eq. (10) has a unique solution. APPENDIX B: OBTAINING INEQUALITIES OF SECTION 2.2 To obtain error estimations for the proposed method, let us rewrite eqs (5), (10) and (11) in the following operator form \begin{eqnarray} {\mathbf {E}}-\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma {\mathbf {E}}={\mathbf {E}}^b, \end{eqnarray} (B1) \begin{eqnarray} \tilde{{\mathbf {E}}}-\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma {\mathbf {U}}={\mathbf {E}}^b, \end{eqnarray} (B2) \begin{eqnarray} {\mathbf {U}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma {\mathbf {U}}=\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}^b. \end{eqnarray} (B3) Note, that operators |$\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma$|, |$\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma$|, |$\left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma \right)^{-1}$|, |$\left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma \right)^{-1}$| are bounded (see Appendices A and C for details). Applying operator |$\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits$| to eq. (B1), we obtain \begin{eqnarray} \mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma {\mathbf {E}}=\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}^b. \end{eqnarray} (B4) Subtracting eq. (B3) from eq. (B1) and adding eq. (B4) gives \begin{eqnarray} {\mathbf {E}}-{\mathbf {U}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma ({\mathbf {E}}-{\mathbf {U}})={\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}. \end{eqnarray} (B5) Therefore, |${\mathbf {E}}-{\mathbf {U}}=\left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma \right)^{-1}({\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}})$| and using eq. (A9), we write \begin{eqnarray} {\left\Vert {{\mathbf {E}}-{\mathbf {U}}}\right\Vert }\le S{\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert }. \end{eqnarray} (B6) Subtracting eq. (B2) from eq. (B1), and using eq. (B5), we obtain \begin{eqnarray} \tilde{{\mathbf {E}}}-{\mathbf {E}}=\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma \left(\mathop {\mathrm{{I}_{}}}\nolimits -\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma \right)^{-1}({\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}). \end{eqnarray} (B7) Hence, \begin{eqnarray} {\left\Vert {\tilde{{\mathbf {E}}}-{\mathbf {E}}}\right\Vert } \le C_E^{\text{in}}{\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert }, \end{eqnarray} (B8) where \begin{eqnarray} C_E^{\text{in}}={\left\Vert {\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma }\right\Vert } S. \end{eqnarray} (B9) The estimation for |${\left\Vert {\tilde{{\mathbf {H}}}-{\mathbf {H}}}\right\Vert }$| can be obtained in a similar manner giving \begin{eqnarray} {\left\Vert {\tilde{{\mathbf {H}}}-{\mathbf {H}}}\right\Vert } \le C_H^{\text{in}}{\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert }, \end{eqnarray} (B10) where \begin{eqnarray} C_H^{\text{in}}={\left\Vert {\mathop {\mathrm{{G}_{H}}}\nolimits \Delta _\sigma }\right\Vert } S. \end{eqnarray} (B11) Estimates for upper bounds for |${\left\Vert {\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma }\right\Vert }$| and |${\left\Vert {\mathop {\mathrm{{G}_{H}}}\nolimits \Delta _\sigma }\right\Vert }$| are presented in Appendix C. Note that the errors in eqs (B8) and (B10) are estimated in |$\mathcal {L}_2$| norm, due to the singularities in integral kernels in eqs (2) and (3) when M ∈ Ω. At the same time for |$M \not\in \mathrm{\Omega }$|, these kernels are smooth, so we can obtain pointwise estimates. To obtain an estimate for |$\left|\tilde{{\mathbf {E}}}(M)-{\mathbf {E}}(M)\right|$| for |$M \not\in \mathrm{\Omega }$|, let us subtract eq. (11) from eq. (2) \begin{eqnarray} {\mathbf {E}}(M)- \tilde{{\mathbf {E}}}(M)=\int \nolimits _{\mathrm{\Omega }}\widehat{G}_E(M,M_0)\Delta _\sigma (M_0)({\mathbf {E}}(M_0)-{\mathbf {U}}(M_0))d\mathrm{\Omega }_{M_0}.\!\!\!\!\!\!\nonumber\\ \end{eqnarray} (B12) Then, for tensor |$\widehat{G}_E(M,M_0)$|, we can write \begin{eqnarray} \widehat{G}_E(M,M_0) \le C_E^{\prime }(M)=\max \limits _{M_0 \in \mathrm{\Omega }}\left|\widehat{G}_E(M,M_0)\right|. \end{eqnarray} (B13) Finally, using the Schwarz inequality, we arrive at the desired estimate \begin{eqnarray} \left|\tilde{{\mathbf {E}}}(M)-{\mathbf {E}}(M)\right| &\le& C_E^{\prime }(M)\max \limits _\mathrm{\Omega }|\Delta _\sigma |\sqrt{V_{\mathrm{\Omega }}}{\left\Vert {{\mathbf {E}}-{\mathbf {U}}}\right\Vert }\nonumber\\ & \le& C_E^{\text{out}}(M){\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert }, \end{eqnarray} (B14) where \begin{eqnarray} C_E^{\text{out}}(M)=C_E^{\prime }(M)\max \limits _\mathrm{\Omega }|\Delta _\sigma |\sqrt{V_{\mathrm{\Omega }}}S, \end{eqnarray} (B15) and where VΩ is the volume of domain Ω. The estimate for |$\left|\tilde{{\mathbf {H}}}(M)-{\mathbf {H}}(M)\right|$| can be obtained similarly \begin{eqnarray} \left|\tilde{{\mathbf {H}}}(M)-{\mathbf {H}}(M)\right|\le C_H^{\text{out}}(M){\left\Vert {{\mathbf {E}}-\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits {\mathbf {E}}}\right\Vert }, \end{eqnarray} (B16) where \begin{eqnarray} C_H^{\text{out}}=C_H^{\prime }(M)\max \limits _\mathrm{\Omega }|\Delta _\sigma |\sqrt{V_{\mathrm{\Omega }}}S \quad C_H^{\prime }(M)=\max \limits _{M_0 \in \mathrm{\Omega }}\left|\widehat{G}_H(M,M_0)\right|.\nonumber\\ \end{eqnarray} (B17) APPENDIX C: MORE ON ESTIMATIONS OF OPERATORS’ NORMS Using eq. (A2), we first write \begin{eqnarray} {\left\Vert {\mathop {\mathrm{{G}_{E}}}\nolimits }\right\Vert } \le \frac{1}{\min \limits _{\mathrm{\Omega }}\mathop {\mathbf {Re}}{{\sigma }}_b}. \end{eqnarray} (C1) Then, taking into account that |${\left\Vert {\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits }\right\Vert }=1$| by design, we have \begin{eqnarray} {\left\Vert {\mathop {\mathrm{{\mathbf {P}}^{N}}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma }\right\Vert } \le {\left\Vert {\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma }\right\Vert } \le \frac{\max \limits _{\mathrm{\Omega }}\left|{\sigma }_a-{\sigma }_b\right|}{\min \limits _{\mathrm{\Omega }}\mathop {\mathbf {Re}}{{\sigma }}_b}. \end{eqnarray} (C2) Obtaining estimates for |${\left\Vert {\mathop {\mathrm{{G}_{H}}}\nolimits }\right\Vert }$| in terms of σa, σb and ω is much more complicated and is beyond the scope of this paper. However, the boundedness can be proved by using the following property of elliptic partial differential equations (cf. Gilbarg & Trudinger 2001) \begin{eqnarray} {\left\Vert {U}\right\Vert }_{W_2^1[\mathrm{\Omega }]}\le C(\mathrm{\Omega },\mathrm{\Omega ^{\prime }}) \left({\left\Vert {U}\right\Vert }_{L_2[\mathrm{\Omega ^{\prime }}]}+{\left\Vert {\mathop {\mathrm{{L}_{}}}\nolimits U}\right\Vert }_{L_2[\mathrm{\Omega ^{\prime }}]}\right), \end{eqnarray} (C3) where |$\mathop {\mathrm{{L}_{ }}}\nolimits$| is the elliptic partial differential operator, Ω ⊃ Ω΄ and C(Ω, Ω΄) > 0 is independent of U. Let us introduce |${\mathbf {V}}=\mathop {\mathrm{{G}_{E}}}\nolimits \Delta _\sigma {\mathbf {U}}$|, then Vγ is the solution of the following equation |$\mathop {\mathrm{{L}_{ }}}\nolimits V_{\gamma }=\Delta _\sigma U_{\gamma }$|, γ = {x, y, z} with corresponding elliptic operator |$\mathop {\mathrm{{L}_{ }}}\nolimits =\mathop {\mathrm{\Delta }}\nolimits +\,i\omega \mu _0{\sigma }_b$| (cf. Ward & Hohmann 1988). Taking into account eqs (C1) and (C3) and expressing |$\mathop {\mathrm{{G}_{H}}}\nolimits$| as |$\mathop {\mathrm{{G}_{H}}}\nolimits =\frac{1}{i\omega \mu _o}\mathop {\mathrm{curl}}\nolimits \mathop {\mathrm{{G}_{E}}}\nolimits$|, we obtain \begin{eqnarray} {\left\Vert {\mathop {\mathrm{{G}_{H}}}\nolimits \Delta _\sigma {\mathbf {U}}}\right\Vert }_{L_2[\mathrm{\Omega }]}&=&\frac{1}{\omega \mu _0}{\left\Vert {\mathop {\mathrm{curl}}\nolimits {\mathbf {V}}}\right\Vert }_{L_2[\mathrm{\Omega }]} \le \frac{1}{\omega \mu _0}{\left\Vert {{\mathbf {V}}}\right\Vert }_{W_2^1[\mathrm{\Omega }]} \nonumber\\ &&\le \frac{1}{\omega \mu _0}C\left({\left\Vert {{\mathbf {V}}}\right\Vert }_{L_2[\mathrm{\Omega ^{\prime }}]}+{\left\Vert {\Delta _\sigma {\mathbf {U}}}\right\Vert }_{L_2[\mathrm{\Omega ^{\prime }}]}\right) \nonumber\\ &&\le \frac{C}{\omega \mu _0}\left({\left\Vert {\mathop {\mathrm{{G}_{E}}}\nolimits }\right\Vert }+1\right)\max \limits _\mathrm{\Omega }\left|{\sigma }_a-{\sigma }_b\right|{\left\Vert {{\mathbf {U}}}\right\Vert }_{L_2[\mathrm{\Omega }]}.\nonumber\\ \end{eqnarray} (C4) Note that we can use Ω in |${\left\Vert {{\mathbf {U}}}\right\Vert }_{L_2[\mathrm{\Omega }]}$| in eq. (C4) instead of Ω΄ because Δσ = 0 outside of Ω. Finally, assuming that |$\mathop {\mathbf {Re}}{{\sigma }}_b(M)>0$|, M ∈ Ω΄ and using eqs (C1) and (C4), we obtain the desired estimate for |${\left\Vert {\mathop {\mathrm{{G}_{H}}}\nolimits \Delta _\sigma }\right\Vert }$| \begin{eqnarray} {\left\Vert {\mathop {\mathrm{{G}_{H}}}\nolimits \Delta _\sigma }\right\Vert }\le C(\mathrm{\Omega },\mathrm{\Omega ^{\prime }})\frac{\max \limits _{\mathrm{\Omega }}\left|{\sigma }_a-{\sigma }_b\right|}{\omega \mu _0}\left(\frac{1}{\min \limits _{\mathrm{\Omega ^{\prime }}}\mathop {\mathbf {Re}}{{\sigma }}_b}+1\right). \end{eqnarray} (C5) © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Using high-order polynomial basis in 3-D EM forward modeling based on volume integral equation method
JF - Geophysical Journal International
DO - 10.1093/gji/ggy059
DA - 2018-05-01
UR - https://www.deepdyve.com/lp/oxford-university-press/using-high-order-polynomial-basis-in-3-d-em-forward-modeling-based-on-fMyYIY1On2
SP - 1387
VL - 213
IS - 2
DP - DeepDyve
ER -