Application of the perfectly matched layer in 3-D marine controlled-source electromagnetic modelling

Application of the perfectly matched layer in 3-D marine controlled-source electromagnetic modelling SUMMARY In this study, the complex frequency-shifted perfectly matched layer (CFS-PML) in stretching Cartesian coordinates is successfully applied to 3-D frequency-domain marine controlled-source electromagnetic (CSEM) field modelling. The Dirichlet boundary, which is usually used within the traditional framework of EM modelling algorithms, assumes that the electric or magnetic field values are zero at the boundaries. This requires the boundaries to be sufficiently far away from the area of interest. To mitigate the boundary artefacts, a large modelling area may be necessary even though cell sizes are allowed to grow toward the boundaries due to the diffusion of the electromagnetic wave propagation. Compared with the conventional Dirichlet boundary, the PML boundary is preferred as the modelling area of interest could be restricted to the target region and only a few absorbing layers surrounding can effectively depress the artificial boundary effect without losing the numerical accuracy. Furthermore, for joint inversion of seismic and marine CSEM data, if we use the PML for CSEM field simulation instead of the conventional Dirichlet, the modelling area for these two different geophysical data collected from the same survey area could be the same, which is convenient for joint inversion grid matching. We apply the CFS-PML boundary to 3-D marine CSEM modelling by using the staggered finite-difference discretization. Numerical test indicates that the modelling algorithm using the CFS-PML also shows good accuracy compared to the Dirichlet. Furthermore, the modelling algorithm using the CFS-PML shows advantages in computational time and memory saving than that using the Dirichlet boundary. For the 3-D example in this study, the memory saving using the PML is nearly 42  per cent and the time saving is around 48 per cent compared to using the Dirichlet. Controlled source electromagnetics (CSEM), Marine electromagnetics, Numerical modelling 1 INTRODUCTION Marine controlled-source electromagnetic (CSEM) method is widely used for investigating the hydrocarbon and gas hydrate resources (Constable 2010). For frequency-domain marine CSEM survey, a towed electric dipole is usually used to generate the EM fields and the recorded data from receivers located at the seafloor can be used to analyse the resistivity distribution of the seafloor. Nowadays 3-D EM modelling studies are of increasing importance for interpreting the EM data acquired in complex geologic settings (Avdeev 2005; Zhdanov 2010). In the 3-D EM modelling problem, the finite-difference/finite-volume, finite-element and integral equation methods are usually used and these methods are well studied in former publications (e.g. Börner 2010; Kelbert et al. 2014). In this study, the staggered finite-difference (SFD) method is used for modelling the 3-D marine CSEM fields, which can produce accurate and reliable numerical solutions with a fast computation rate (Mackie et al. 1993; Smith 1996; Weiss & Constable 2006; Sasaki & Meju 2009; Mittet 2010). For these traditional frameworks of 3-D EM modelling algorithms, a Dirichlet boundary is usually used, which assumes the electric or magnetic field values are zero at the boundaries. This is a truncated boundary and it requires the boundaries to be sufficiently far away from the area of interest. To mitigate the boundary artefacts, a large modelling area may be necessary though cell sizes are allowed to grow toward the boundaries due to the diffusion of the electromagnetic wave propagation. An alternative way is to use the perfectly matched layer (PML), which is the common choice for seismic wave propagating simulation (e.g. Hu et al. 2007; Pan et al. 2012). By using the PML, the modelling area of interest could be restricted to the target region and only a few surrounding absorbing layers can effectively depress the artificial boundary effect without losing the numerical accuracy. It is worth mentioning that, for joint inversion of marine CSEM and seismic data, if we use the PML for CSEM field simulation instead of the conventional Dirichlet, the modelling area for these two different geophysical data collected from the same survey area could be the same, which is convenient for joint inversion grid matching (Hu et al. 2009). Chen et al. (1997) applied Bérenger’s PML to 2-D transient EM modelling problem using the finite-difference time-domain (FDTD) modelling of lossless media (Bérenger 1994). Schwarzbach et al. (2011) successfully utilized Bérenger’s PML in 3-D frequency-domain CSEM diffusing modelling problem with adaptive higher order finite element method. Similarly, Mittet (2010), de la Kethulle de Ryhove & Mittet (2014) used Bérenger’s PML in 3-D marine magnetotelluric modelling by a finite-difference time-domain algorithm in which Maxwell’s equations are solved in a fictitious-wave domain. The original PML proposed by Bérenger (1994) can be optimized and simply implemented in the form of the complex coordinate stretching along Cartesian coordinates (Chew & Weedon 1994; Chew & Jin 1996). To deduce the asymptotic behaviours at low and high frequencies, a complex frequency-shifted PML (CFS-PML) is also developed which greatly improve the absorbing performance of PML (Kuzuoglu & Mittra 1996; Bérenger 2002). Hu et al. (2017) applied the CFS-PML boundary condition for transient electromagnetic modelling using a fictitious wave domain method. In this study, we present a novel 3-D frequency-domain marine CSEM modelling algorithm, in which CFS-PML is successfully applied. This study can be divided into three parts. We first give an introduction of 3-D frequency-domain CSEM modelling scheme using the SFD grids with the CSF-PML. To avoid the source singularities, the secondary-field approach is used and the primary fields excited by the electric dipole source could be calculated quasi-analytically for the 1-D layered background. Then we compare the performance of CFS-PML with the conventional Dirichlet boundary by numerical analysis. Finally we give the conclusion that the CFS-PML modelling scheme proposed shows advantages in computational time and memory saving compared to using the Dirichlet boundary without losing numerical accuracy. 2 FORMULATIONS 2.1 Governing equations The electric/magnetic fields can be split into a primary part and a secondary (or scattered) part to avoid source singularities (Newman & Alumbaugh 1995; Sasaki & Meju 2009; Streich 2009). The primary fields are excited by the arbitrarily oriented dipole sources and can be computed quasi-analytically following Li & Li (2016). The secondary fields are computed by the frequency-domain SFD method (Yee 1966). Assuming that the time factor is e − iωt where i2 = −1, by using the electric dipole source, the Maxwell’s equations in the frequency domain can be written as   \begin{eqnarray} \nabla \times {{\bf E}}^S - {\rm i}\omega \mu _0 {{\bf H}}^S = 0, \quad \nabla \times {{\bf H}}^S - \sigma ^\ast {{\bf E}}^S = (\sigma ^\ast - \sigma ^{P\ast }) {{\bf E}}^P, \nonumber\\ \end{eqnarray} (1)where superscripts P and S are for the primary and secondary fields, respectively. E and H are the electric and magnetic fields, ω is the angular frequency, μ0 is the magnetic permeability in free space, and the complex conductivity σ* = σ − iωε0 consists of the electric conductivity σ and permittivity ε0 in free space. The secondary-source term is given by (σ* − σP*)EP where σ* = σP* + σS*. After eliminating HS, eq. (1) becomes   \begin{eqnarray} \nabla \times \nabla \times {{\bf E}}^S - {\rm i}\omega \mu _0 \sigma ^\ast {{\bf E}}^S = {\rm i}\omega \mu _0 (\sigma ^\ast - \sigma ^{P\ast }) {{\bf E}}^P. \end{eqnarray} (2)In this study, we use eq. (2) for simulating marine CSEM fields with the frequency range of 0.1–10 Hz. 2.2 Implementation of CFS-PML In the complex stretched coordinate approach (Chew & Weedon 1994; Kuzuoglu & Mittra 1996), the governing equation using CFS-PML can be rewritten as   \begin{eqnarray} \nabla _h \times \nabla _e \times {{\bf E}}^S - {\rm i}\omega \mu _0 \sigma ^\ast {{\bf E}}^S = {\rm i}\omega \mu _0 (\sigma ^\ast - \sigma ^{P\ast }) {{\bf E}}^P. \end{eqnarray} (3)where $$\nabla _e = (\frac{1}{\gamma ^e_x} \partial _x,\frac{1}{\gamma ^e_y} \partial _y,\frac{1}{\gamma ^e_z} \partial _z)$$, $$\nabla _h = (\frac{1}{\gamma ^h_x} \partial _x,\frac{1}{\gamma ^h_y} \partial _y,\frac{1}{\gamma ^h_z} \partial _z)$$, $$\gamma ^e_\nu$$ and $$\gamma ^h_\nu$$ (ν = x, y, z) are the complex stretched coordinate factors (decay factors) for electric and magnetic field components, respectively. SFD discretization of eq. (3) is similar to that in Newman & Alumbaugh (1995) and Streich et al. (2013) (Appendix A). The entire computational domain is divided into Ncell = Nx × Ny × Nz cells, where Nx, Ny, and Nz are the number of cells in the x-, y-, and z- directions, respectively. After discretization, eq. (3) can be written in the matrix form as   \begin{equation} {{\bf K}} {{\bf U}} = {{\bf P}}, \end{equation} (4)where K is a symmetric complex matrix of dimension (3Nx × Ny × Nz)2; the unknown vector U of length 3Nx × Ny × Nz contains the electric field values $$E_x^S$$, $$E_y^S$$ and $$E_z^S$$ for all nodes; and P of length 3Nx × Ny × Nz is the secondary source vector of the right-hand side of eq. (4) given by eq. (3). The entries in K depend on the grid spacing and the frequency-dependent medium properties, and K is up to 13 nonzero entries per line. The 3-D array electric field values $$E_{i+\frac{1}{2},j,k}^{xS}$$, $$E_{i,j+\frac{1}{2},k}^{yS}$$, and $$E_{i,j,k+\frac{1}{2}}^{zS}$$, where the subscripts i = 1, …, Nx, j = 1, …, Ny, and k = 1, …, Nz, can be mapped into a 1-D column array of indices (i + 1/2, j, k) → [(i + 1/2 − 1) × Ny + j] × Nz + k, (i, j + 1/2, k) → [(i − 1) × Ny + j + 1/2] × Nz + k, and (i, j, k + 1/2) → [(i − 1) × Ny + j] × Nz + k + 1/2, respectively. Similarly, we can write the source term in a 1-D column array. The stiffness matrix K is symmetric and highly sparse. Fig. 1 shows an example of the symmetric structure of the stiffness matrix K where a 5 × 5 × 5 staggered grid is used. Figure 1. View largeDownload slide Illustration of the symmetric structure of the stiffness matrix for a 5 × 5 × 5 staggered grid. Note that only the locations of nonzero entries before applying the Dirichlet or PML boundary condition are shown (blue points). nz = 1944 is the number of nonzero entries of the stiffness matrix. Figure 1. View largeDownload slide Illustration of the symmetric structure of the stiffness matrix for a 5 × 5 × 5 staggered grid. Note that only the locations of nonzero entries before applying the Dirichlet or PML boundary condition are shown (blue points). nz = 1944 is the number of nonzero entries of the stiffness matrix. The formed linear equations in (4) given by the discretization of eq. (3) are solved by a multifrontal direct solver MUMPS 5.0.2 parallelized by OpenMP (Amestoy et al. 2001, 2012), which could avoid uncertainties in pre-conditioning and convergence for iterative solutions, especially for low frequencies (Farquharson & Miensopust 2011; Oldenburg et al. 2013). In this section, we will focus on the implementation of CFS-PML in detail. Following Chew & Weedon (1994) and Bérenger (2007b), the PML decay factors are given as   \begin{eqnarray} \gamma ^e_\nu &=& 1 - {\rm i} \eta ^e_\nu / \omega , \quad \nu =x,y,z \end{eqnarray} (5)  \begin{eqnarray} \gamma ^h_\nu &=& 1 - {\rm i} \eta ^h_\nu / \omega , \quad \nu =x,y,z. \end{eqnarray} (6)Absorbing boundaries at the edges of the simulation region may be created by choosing appropriate values of ην. Usually, the decay factor ην varies gradually from 0 at the PML/non-PML interface to its maximum at the outer boundaries of the computational domain to minimize the numerical reflections caused by spatial discretization (Hu et al. 2007). In this study, the PML decay factor ην is determined empirically using the polynomial form of the PML (Hu et al. 2007). For example, the artificial attenuation along the x-direction is defined as (Bérenger 1996; Pan et al. 2012):   \begin{eqnarray} \eta _x = \left\lbrace \begin{array}{rcl}\eta _{\rm {max}} \left( -\frac{x}{L_{\rm {pml}}} \right)^m, & &\quad {{\rm if}\,\quad -L_{\rm {pml}} \le x < 0}\\ {0,} & &\quad {{\rm if}\,\quad 0 \le x \le 1}\\ {\eta _{\rm {max}} \left( \frac{x-1}{L_{\rm {pml}}} \right)^m,} & &\quad {{\rm if}\,\quad 1 < x \le 1+L_{\rm {pml}}} \end{array} \right. \end{eqnarray} (7)where m is either 2 or 3, and this value depends on the size of the computational domain. As pointed out by Bérenger (2002), m = 3 has a better absorption effect on the evanescent region than m = 2, whereas m = 2 has a better absorption effect on the traveling region. We use m = 2 in this study. The length of the interior domain is scaled to 1, Lpml is the thickness of the PML layer on each side, ηmax is the maximum PML decay factor. We use the same definitions for the artificial attenuation parameters in the y- and z-directions. Fig. 2 shows the different PML layers surrounding the original computational area. Note for the computation of the solution inside the PML layers in the x-direction, ηx is positive and ηy, ηz = 0. For PML layers in the y-direction, ηy is positive and ηx, ηz = 0. For PML layers in the z-direction, ηz is positive and both ηx, ηy = 0 (see Table 1). For PML layers in the corners, where we need damping in all directions, ηx, ηy and ηz are all positive (Collino & Tsogka 2001). Figure 2. View largeDownload slide The schematic map showing different PML layers surrounding the 3-D original computational cube. The left plot shows the original computational volume, which is a cuboid without PML layers. The right plot shows the PML layers surrounding the surrounding the original computational volume. The PML layers can be divided into seven different types, which have different values of ηx, ηy and ηz (see Table 1). Figure 2. View largeDownload slide The schematic map showing different PML layers surrounding the 3-D original computational cube. The left plot shows the original computational volume, which is a cuboid without PML layers. The right plot shows the PML layers surrounding the surrounding the original computational volume. The PML layers can be divided into seven different types, which have different values of ηx, ηy and ηz (see Table 1). Table 1. The damping factors ην(ν = x, y, z) for different PML layers. PML layers  PML-x  PML-y  PML-z  PML-x, y  PML-y, z  PML-z, x  PML-x, y, z  ηx  >0  =0  =0  >0  =0  >0  >0  ηy  =0  >0  =0  >0  >0  =0  >0  ηz  =0  =0  >0  =0  >0  >0  >0  PML layers  PML-x  PML-y  PML-z  PML-x, y  PML-y, z  PML-z, x  PML-x, y, z  ηx  >0  =0  =0  >0  =0  >0  >0  ηy  =0  >0  =0  >0  >0  =0  >0  ηz  =0  =0  >0  =0  >0  >0  >0  View Large The complex frequency stretched (CFS) PML, as introduced in Kuzuoglu & Mittra (1996), Roden & Gedney (2000) and Pan et al. (2012), is one of a number of approaches that are suitable to mimic the reflection-free boundaries at the outermost surface of the interior domain. The complex frequency stretched coordinate factors are given by   \begin{eqnarray} \gamma ^e_\nu = \kappa ^e_\nu + \frac{\eta ^e_\nu }{\alpha ^e_\nu + \rm {i}\omega }, \end{eqnarray} (8)  \begin{eqnarray} \gamma ^h_\nu = \kappa ^h_\nu + \frac{\eta ^h_\nu }{\alpha ^h_\nu + \rm {i}\omega }, \end{eqnarray} (9)where ην and κν are positive real and κν is ≥1 (Bérenger 2007a). ην are the PML decay factors, αν are point-wise coefficients related to the CFS-PML (Roden & Gedney 2000). The essence of the real parameter αν is to absorb the evanescent waves and to improve the absorption performance at grazing angles (Bérenger 2002; Komatitsch & Martin 2007; Martin et al. 2009). In this study, the parameter κν is set to be 1, since the effect of these parameters for the absorption is not significant (Bérenger 2002). In the PML region, the conductivity value is equal to that of the outmost cell of the interior domain. This strategy is employed to minimize the numerical reflection error at the interface between the PML domain and the interior domain. We use similar setting for αν as Roden & Gedney (2000). We choose to make αν vary in a linear fashion in their respective PML layer between a maximum value αmax at the beginning (i.e. the entrance) of the PML and zero at its top. The maximum value of αν is set to be 95 per cent of the transmitting frequency. With these settings, the discretized eq. (3) has the same formulas everywhere in space. The only difference between the PML domain and the interior domain equation is the value of the artificial attenuation ην. To reduce the discretization error, the parameter ην are scaled such that they are 0 and 1 at the PML/interior volume interface, respectively, and are maximum at the exterior boundary. We use an optimal maximum attenuation (ηmax) formulation similar to that in Komatitsch & Martin (2007), Martin et al. (2009) and Pan et al. (2012)   \begin{equation} \eta _{\rm {max}} (x,y,z) = c \frac{(m+1)\sigma (x,y,z)}{L_{\rm {pml}}}, \end{equation} (10)where σ(x, y, z) is the conductivity of the cell adjacent to the PML grid, Lpml is the thickness of the PML domain and c a constant that determines the reflection from the PML interface. In this study, the empirical parameter c is chosen to be several times of the skin depth for the background model. 3 NUMERICAL ANALYSIS In this section, we examine the performance of the CFS-PML proposed. Two examples are tested, that are, the 1-D case compared with the quasi-analytical solutions given by Li & Li (2016) and the 3-D case compared with the solutions of adaptive finite-element method (adaptive FEM) given by MARE3DEM code (Zhang 2017). We have written a Fortran 90 code to implement the 3-D SFD modelling scheme for marine CSEM data. The numerical accuracy, time consumption and memory saving are compared between the CFS-PML and the conventional Dirichlet boundary (Li et al. 2017). The numerical test is performed under the Dell Precision Tower 3620 (3.50 GHz CPU Intel@Xeon E3-1240 v5 family including two processors with up to four cores per processor, memory up to 16 G), which is suitable for OpenMP parallel programming. Note the MUMPS 5.0.2 used for factorizing the complex sparse matrix formed by SFD discretization is parallelized by OpenMP in an ‘out-of-core’ environment (Amestoy et al. 2001, 2012). When the ‘out-of-core’ phase is activated and the complete matrix of factors is written to disk and will be read each time a solution phase is requested, therefore the memory requirement can be significantly reduced while not increasing much of the factorization time on a reasonably small number of processors. 3.1 The 1-D case For simplicity, we first use a 1-D canonical reservoir model given by Constable & Weiss (2006), so that our SFD numerical solutions can be easily validated by the results obtained from the quasi-analytical 1-D algorithm (Li & Li 2016). The 1-D canonical reservoir model consists of a sea-water layer with 0.3 Ωm resistivity and 1 km thickness, a sediment layer with 1 Ωm resistivity and a 1000 m thickness (see Fig. 3). The canonical reservoir is buried at a depth of 2 km below the sea surface with 100 Ωm resistivity and a 100 m thickness. Inline (transmitter pointing along y-direction) transmissions are excited at a frequency of 0.25 Hz. The transmitter location is x = 0, y = −5000, z = 950 m, and 51 receivers are located at the seafloor with 200-m intervals from y = −5000 to 5000 m along the towed line. Here, we only consider the electric and magnetic fields in the transmitter-receiver plane at x = 0. Figure 3. View largeDownload slide A 1-D canonical reservoir model similar to Constable & Weiss (2006). There are five layers including air layer, sea-water layer and three seafloor layers. The • and ▾ indicate the transmitter and receivers, respectively. Figure 3. View largeDownload slide A 1-D canonical reservoir model similar to Constable & Weiss (2006). There are five layers including air layer, sea-water layer and three seafloor layers. The • and ▾ indicate the transmitter and receivers, respectively. The computational volume is {(x, y, z): −50 km ≤ x, y, z ≤ 50 km} and it is divided into 76 × 76 × 72 staggered-grids when using the Dirichlet boundary condition (Figs 4a). The x- and y- grid spacings are the same and becomes larger towards the boundaries, and the z-grid spacing becomes larger with depth. The minimum horizontal spacing is set to be 200 m, while the minimum vertical grid spacing is 50 m. For the CFS-PML, the computational volume is {(x, y, z): −6.6 km ≤ x, y ≤ 6.6 km, −8 km ≤ z ≤ 3 km} and a 66 × 66 × 64 staggered-grid including 8 PML layers surrounding the computational area is used (Fig. 4b). Figure 4. View largeDownload slide A 2-D cutaway of the mesh in the transmitter–receiver plane at x = 0 used for forward calculation of the 1-D canonical reservoir model in Fig. 3 using a 3-D SFD code with a PML and Dirichlet boundary condition. (a) The mesh using the conventional Dirichlet boundary in which local mesh for reservoir area bounded by yellow box is enlarged. (b) Local mesh using the CFS-PML boundary. Figure 4. View largeDownload slide A 2-D cutaway of the mesh in the transmitter–receiver plane at x = 0 used for forward calculation of the 1-D canonical reservoir model in Fig. 3 using a 3-D SFD code with a PML and Dirichlet boundary condition. (a) The mesh using the conventional Dirichlet boundary in which local mesh for reservoir area bounded by yellow box is enlarged. (b) Local mesh using the CFS-PML boundary. Fig. 5 shows the 3-D SFD solutions compared to 1-D quasi-analytic solutions for the frequency 0.25 Hz. From Fig. 5 we can see that, the numerical accuracy is acceptable for both the gridding using Dirichlet and CFS-PML boundary. The absolute errors of the amplitude for Ey, Ez and Hx are no more than 2 per cent, 4 per cent and 2 per cent, respectively. The relative errors of the phase for Ey, Ez and Hx are all less than 1°. When using the CFS-PML boundary, less grids are used so that both the memory and computation time saving are significant (Table 2). For this 1-D test, the memory saving is nearly 43  per cent and the time saving is around 50 per cent. Figure 5. View largeDownload slide 3-D numerical solutions against 1-D quasi-analytic solutions for the frequency 0.25 Hz. The y-axis in km is also the source receiver offset. The SFD solutions with the CFS-PML boundary are compared to those with the Dirichlet boundary condition given by Li et al. (2017). (a) Amplitude of the electric field E and electric field H in which the units are VA − 1 m − 2 and m − 2, respectively. (b) Phase in degrees for the electric field E and electric field H. (c) Absolute error of amplitude for the electric field E and electric field H compared to the 1-D analytic solution given by Li & Li (2016). (d) Relative error of phase for the electric field E and electric field H compared to the 1-D analytic solution. The solid line indicates the quasi-analytic results for the 1-D model in Fig. 3, × and ◯ indicate the SFD results using the Dirichlet boundary and the PML boundary, respectively. The colours black, blue and red are for Ey, Ez and Hx, respectively. Figure 5. View largeDownload slide 3-D numerical solutions against 1-D quasi-analytic solutions for the frequency 0.25 Hz. The y-axis in km is also the source receiver offset. The SFD solutions with the CFS-PML boundary are compared to those with the Dirichlet boundary condition given by Li et al. (2017). (a) Amplitude of the electric field E and electric field H in which the units are VA − 1 m − 2 and m − 2, respectively. (b) Phase in degrees for the electric field E and electric field H. (c) Absolute error of amplitude for the electric field E and electric field H compared to the 1-D analytic solution given by Li & Li (2016). (d) Relative error of phase for the electric field E and electric field H compared to the 1-D analytic solution. The solid line indicates the quasi-analytic results for the 1-D model in Fig. 3, × and ◯ indicate the SFD results using the Dirichlet boundary and the PML boundary, respectively. The colours black, blue and red are for Ey, Ez and Hx, respectively. Table 2. Comparison between the CFS-PML boundary and the conventional Dirichlet boundary for the 1-D example. Boundary  Gridding  PML layers  Number of unknowns  Memory used (G)  Total time used (min)  Dirichlet  76 × 76 × 72  –  1247616  ≈7.76  ≈14.15  CFS-PML  66 × 66 × 64  8  836352  ≈4.41  ≈7.04  Boundary  Gridding  PML layers  Number of unknowns  Memory used (G)  Total time used (min)  Dirichlet  76 × 76 × 72  –  1247616  ≈7.76  ≈14.15  CFS-PML  66 × 66 × 64  8  836352  ≈4.41  ≈7.04  View Large Similar to Roden & Gedney (2000) and Hu et al. (2017), the reflection errors at the receiver positions are calculated using   \begin{eqnarray} {\rm Reflection\,\,error} = 20\log _{10} \left| \frac{{\rm abs}({{\bf F}}_i) - {\rm abs}({{\bf F}}^0_i)}{{\rm abs}({{\bf F}}^0_i)} \right|, i=1,\ldots ,N_{R_x} \nonumber\\ \end{eqnarray} (11)where $$N_{R_x}$$ is the number of receivers, F represents the electric E or magnetic field H observed at the receiver locations using the CFS-PML boundary, F0 represents the reference field with boundaries distant enough to avoid any boundary reflections that might interfere with the observed data. In this example, the 1-D quasi-analytical results are taken as the reference fields. Fig. 6 shows the reflection error for both the electric field (Fig. 6a) and the magnetic field (Fig. 6b) at the receiver positions for the 1-D model shown in Fig. 3. Actually, the reflection errors at the receiver positions in eq. (11) indicate the total electric or magnetic field reflection level. The reflection error of both E and H using CFS-PML boundary is low (less than −30 dB noise which corresponds to the relative error of 3 per cent), which is similar to that using the Dirichlet boundary. This indicates that the CFS-PML boundary could effectively depress the artificial boundary effect. Figure 6. View largeDownload slide Reflection error of the PML and Dirichlet boundaries for the 1-D model shown in Fig. 3. The symbols × and ◯ are for the Dirichlet boundary and the PML boundary, respectively. (a) Reflection error for the electric field E, and (b) for the magnetic field H. Figure 6. View largeDownload slide Reflection error of the PML and Dirichlet boundaries for the 1-D model shown in Fig. 3. The symbols × and ◯ are for the Dirichlet boundary and the PML boundary, respectively. (a) Reflection error for the electric field E, and (b) for the magnetic field H. 3.2 The 3-D case In this part, we present a 3-D example similar to that given by Weiss & Constable (2006) (Fig. 7). The 3-D geo-electric model consists of a reservoir represented by a 100-m thick horizontal slab and resistivity 100 Ωm. The reservoir body is buried at a depth of 1 km within a half-space of 1 Ωm water-saturated sediment. The thin slab is a cuboid with 4 km length, 4 km width and 100 m height. The centre of the top surface of the slab corresponds to the origin o of the Cartesian coordinates used. Inline (transmitter pointing along y-direction) transmissions are excited at a frequency of 1 Hz. The transmitter location is x, y = 0, z = 900 m, and 21 receivers are located at the seafloor with 200 m intervals from y = 0 to 5000 m along the towed line. Only the electric and magnetic fields in the transmitter-receiver plane at x = 0 are considered. Figure 7. View largeDownload slide A 3-D canonical reservoir model similar to Weiss & Constable (2006). There are air layer, sea-water layer and a canonical slab of the volume {(x, y, z): −2 km ≤ x, y ≤ 2 km, 2 km ≤ z ≤ 2.1 km} which is embedded into a seafloor half-space. The • and ▾ indicate the transmitter and receivers, respectively. Figure 7. View largeDownload slide A 3-D canonical reservoir model similar to Weiss & Constable (2006). There are air layer, sea-water layer and a canonical slab of the volume {(x, y, z): −2 km ≤ x, y ≤ 2 km, 2 km ≤ z ≤ 2.1 km} which is embedded into a seafloor half-space. The • and ▾ indicate the transmitter and receivers, respectively. The computational volume is {(x, y, z): −50 km ≤ x, y, z ≤ 50 km} is also divided into 78 × 78 × 74 staggered-grids when using the Dirichlet boundary condition (Figs 8a). For the CFS-PML, the computational volume is {(x, y, z): −6.6 ≤ x, y ≤ 6.6 km, −8 km ≤ z ≤ 3 km} and a 68 × 68 × 64 staggered-grid including 8 PML layers surrounding the computational area is used (Fig. 8b). For comparison, the adaptive FEM solutions computed by MARE3DEM are used (Zhang 2017). The computational volume is {(x, y, z): −50 ≤ x, y, z ≤ 50 km} and the total triangular elements are 755 804 after 10 refinement iterations. Figure 8. View largeDownload slide A 2-D cutaway of the mesh in the transmitter-receiver plane at x = 0 used for forward calculation of the 3-D canonical reservoir model in Fig. 7 using a 3-D SFD code with a PML and Dirichlet boundary condition. (a) The mesh using the conventional Dirichlet boundary, in which local mesh for reservoir area bounded by yellow box is enlarged. (b) Local mesh using the CFS-PML boundary. Figure 8. View largeDownload slide A 2-D cutaway of the mesh in the transmitter-receiver plane at x = 0 used for forward calculation of the 3-D canonical reservoir model in Fig. 7 using a 3-D SFD code with a PML and Dirichlet boundary condition. (a) The mesh using the conventional Dirichlet boundary, in which local mesh for reservoir area bounded by yellow box is enlarged. (b) Local mesh using the CFS-PML boundary. Fig. 9 shows the 3-D SFD solutions compared to 3-D adaptive FEM solutions for the frequency 0.25 Hz. Fig. 9 shows a good numerical accuracy for using Dirichlet or CFS-PML boundary. The relative errors of the amplitude for Ey, Ez and Hx are all no more than 1.2 per cent. The absolute errors of the phase are all less than 1°. Both the memory and computation time saving are significant when using the CFS-PML boundary with less grids (Table 3). For the 3-D test, the memory saving using the PML is nearly 42  per cent and the time saving is around 48 per cent compared to using the Dirichlet. In Fig. 9, one can distinguish the 3-D adaptive finite-element results (solid line) and quasi-analytic results for the 1-D background model without the 3-D resistive slab (dashed line). This indicates the effect of the 3-D thin slab is apparent. Figure 9. View largeDownload slide 3-D numerical solutions against 3-D adaptive FEM solutions of MARE3DEM given by Zhang (2017) for the frequency 1 Hz. The SFD solutions with the CFS-PML boundary are compared to those with the Dirichlet boundary condition given by Li et al. (2017). Note that the thin slab is of dimension −4 km to 4 km along y axis. The solid line indicates the adaptive finite-element results, the dash line indicates the quasi-analytic results for the 1-D background model in Fig. 7 without the 3-D resistive slab, × and ◯ indicate the SFD results using the Dirichlet boundary and the PML boundary, respectively. The colours black, blue and red are for Ey, Ez and Hx, respectively. Figure 9. View largeDownload slide 3-D numerical solutions against 3-D adaptive FEM solutions of MARE3DEM given by Zhang (2017) for the frequency 1 Hz. The SFD solutions with the CFS-PML boundary are compared to those with the Dirichlet boundary condition given by Li et al. (2017). Note that the thin slab is of dimension −4 km to 4 km along y axis. The solid line indicates the adaptive finite-element results, the dash line indicates the quasi-analytic results for the 1-D background model in Fig. 7 without the 3-D resistive slab, × and ◯ indicate the SFD results using the Dirichlet boundary and the PML boundary, respectively. The colours black, blue and red are for Ey, Ez and Hx, respectively. Table 3. Comparison between the CFS-PML boundary and the conventional Dirichlet boundary for the 3-D example. Boundary  Gridding  PML layers  Number of unknowns  Memory used (G)  Total time used (min)  Dirichlet  78 × 78 × 74  –  1350648  ≈8.75  ≈15.46  CFS-PML  68 × 68 × 64  8  887808  ≈5.09  ≈8.08  Boundary  Gridding  PML layers  Number of unknowns  Memory used (G)  Total time used (min)  Dirichlet  78 × 78 × 74  –  1350648  ≈8.75  ≈15.46  CFS-PML  68 × 68 × 64  8  887808  ≈5.09  ≈8.08  View Large We notice that the PML appears to have a little higher error compared to the Dirichlet for both the 1-D and 3-D examples. We think that because the modelling area is large enough and the boundaries are far away enough (the modelling area is set to be {(x, y, z): −50 km ≤ x, y, z ≤ 50 km}), no or little boundary reflections appear for the Dirichlet. For the PML, the modelling area for PML is still much smaller than that of Dirichlet. Although we use 10 PML layers which are proved to be sufficient enough for absorbing boundary reflections, the absorbing rate cannot reach 100 per cent. The reflection errors at the receiver positions for the 3-D example are shown in Fig. 10. The 3-D adaptive FEM results, using the Dirichlet boundary, are taken as the reference fields. In this test, the Dirichlet boundary for finite element simulation is far enough to avoid any boundary reflections that might interfere with the observed data. The reflection errors of both E and H using CFS-PML boundary are less than −30 dB and −28 dB, respectively, which are similar to those using the Dirichlet boundary. This indicates that the artificial boundary effect is effectively depressed by the CFS-PML boundary. Figure 10. View largeDownload slide Reflection error of the PML and Dirichlet boundaries for the 3-D model shown in Fig. 7. Note that the thin slab is of dimension −4 km to 4 km along y axis. The symbols × and ◯ are for the Dirichlet boundary and the PML boundary, respectively. (a) Reflection error for the electric field E, and (b) for the magnetic field H. Figure 10. View largeDownload slide Reflection error of the PML and Dirichlet boundaries for the 3-D model shown in Fig. 7. Note that the thin slab is of dimension −4 km to 4 km along y axis. The symbols × and ◯ are for the Dirichlet boundary and the PML boundary, respectively. (a) Reflection error for the electric field E, and (b) for the magnetic field H. 4 CONCLUSIONS In this study, the CFS-PML has been applied to the 3-D marine CSEM forward problem. For the PML boundary, the model area can be restricted to the region of interest and only a few absorbing layers surrounding can effectively depress the artificial boundary effect without losing the numerical accuracy. The next step is to develop an inversion scheme for interpreting real data based on the 3-D marine CSEM modelling algorithm using CFS-PML. Furthermore, for developing a joint 3-D inversion scheme of marine CSEM and seismic data, the proposed 3-D CSEM modelling scheme using CFS-PML could be more convenient than using the conventional Dirichlet boundary condition. For seismic wave simulation, a PML boundary is commonly used. If we use the PML for CSEM field simulation, the modelling area for these two different geophysical data collected from the same survey area could be the same. This could avoid dealing with the grid matching when the CSEM gridding area using Dirichlet boundary is usually larger than the seismic gridding area (Hu et al. 2009). ACKNOWLEDGEMENTS The authors acknowledge the editor Jörg Renner, the associate editor Mark Everett, and two anonymous reviewers for their constructive comments. Wenyi Hu is acknowledged for his fruitful discussions about the implementation of PML. Yuxiang Zhang is thanked for providing the 3-D adaptive finite-element results using MARE3DEM. Gang Li thanks the support from the Sino-German (CSC-DAAD) Postdoc Scholarship Program 2016. This work is funded by the National Natural Science Foundation of China (41774080, 41704075 and 41604063) and the Shandong Provincial Natural Science Foundation, China (ZR2016DQ15 and ZR2016DB31). Finally, Gang Li appreciates the support and help from Marion Jegen and Sebastian Hoelz for his study at GEOMAR. REFERENCES 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   Amestoy P.R., Duff I.S., L’Excellent J.-Y., Robert Y., Rouet F.-H., Uçar B., 2012. On computing inverse entries of a sparse matrix in an out-of-core environment, SIAM J. Sci. Comput. , 34( 4), A1975– A1999. Google Scholar CrossRef Search ADS   Avdeev D.B., 2005. Three-dimensional electromagnetic modelling and inversion from theory to application, Surv. Geophys. , 26( 6), 767– 799. Google Scholar CrossRef Search ADS   Bérenger J.-P., 1994. A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. , 114( 2), 185– 200. Google Scholar CrossRef Search ADS   Bérenger J.-P., 1996. Three-dimensional perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. , 127( 2), 363– 379. Google Scholar CrossRef Search ADS   Bérenger J.-P., 2002. Application of the CFS PML to the absorption of evanescent waves in waveguides, IEEE Microw. Wirel. Compon. Lett. , 12( 6), 218– 220. Google Scholar CrossRef Search ADS   Bérenger J.-P., 2007a. On the Huygens absorbing boundary conditions for electromagnetics, J. Comput. Phys. , 226( 1), 354– 378. Google Scholar CrossRef Search ADS   Bérenger J.-P., 2007b. Perfectly Matched Layer (PML) for Computational Electromagnetics , 1st edn, Morgan & Claypool. Börner R.-U., 2010. Numerical modelling in geo-electromagnetics: advances and challenges, Surv. Geophys. , 31( 2), 225– 245. Google Scholar CrossRef Search ADS   Chen Y.H., Chew W.C., Oristaglio M.L., 1997. Application of perfectly matched layers to the transient modeling of subsurface EM problems, Geophysics , 62( 6), 1730– 1736. Google Scholar CrossRef Search ADS   Chew W.C., Jin J.M., 1996. Perfectly matched layers in the discretized space: an analysis and optimization, Electromagnetics , 16( 4), 325– 340. Google Scholar CrossRef Search ADS   Chew W.C., Weedon W.H., 1994. A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates, Microw. Opt. Technol. Lett. , 7( 13), 599– 604. Google Scholar CrossRef Search ADS   Collino F., Tsogka C., 2001. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media, Geophysics , 66( 1), 294– 307. 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., Weiss C.J., 2006. Mapping thin resistors and hydrocarbons with marine EM methods: insights from 1D modeling, Geophysics , 71( 2), G43– G51. Google Scholar CrossRef Search ADS   de la Kethulle de Ryhove S., Mittet R., 2014. 3D marine magnetotelluric modeling and inversion with the finite-difference time-domain method, Geophysics , 79( 6), E269– E286. Google Scholar CrossRef Search ADS   Farquharson C.G., Miensopust M.P., 2011. Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction, J. Appl. Geophys. , 75( 4), 699– 710. Google Scholar CrossRef Search ADS   Haber E., Ascher U., Aruliah D., Oldenburg D., 2000. Fast simulation of 3D electromagnetic problems using potentials, J. Comput. Phys. , 163( 1), 150– 171. Google Scholar CrossRef Search ADS   Hu W., Abubakar A., Habashy T.M., 2007. Application of the nearly perfectly matched layer in acoustic wave modeling, Geophysics , 72( 5), SM169– SM175. Google Scholar CrossRef Search ADS   Hu W., Abubakar A., Habashy T.M., 2009. Joint electromagnetic and seismic inversion using structural constraints, Geophysics , 74( 6), R99– R109. Google Scholar CrossRef Search ADS   Hu Y., Egbert G., Ji Y., Fang G., 2017. A novel CFS-PML boundary condition for transient electromagnetic simulation using a fictitious wave domain method, Radio Sci. , 52( 1), 118– 131. Google Scholar CrossRef Search ADS   Kelbert A., Kuvshinov A., Velimsky J., Koyama T., Ribaudo J., Sun J., Martinec Z., Weiss C.J., 2014. Global 3-D electromagnetic forward modelling: a benchmark study, Geophys. J. Int. , 197( 2), 785– 814. Google Scholar CrossRef Search ADS   Komatitsch D., Martin R., 2007. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation, Geophysics , 72( 5), SM155– SM167. Google Scholar CrossRef Search ADS   Kuzuoglu M., Mittra R., 1996. Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers, IEEE Microw. Guid. Wave Lett. , 6( 12), 447– 449. Google Scholar CrossRef Search ADS   Li Y., Li G., 2016. Electromagnetic field expressions in the wavenumber domain from both the horizontal and vertical electric dipoles, J. Geophys. Eng. , 13( 4), 505– 515. Google Scholar CrossRef Search ADS   Li G., Zhang L., Han B., 2016. Stable electromagnetic modeling using a multigrid solver on stretching grids: The magnetotelluric example, IEEE Geosci. Remote Sens. Lett. , 13, 334– 338. Li G., Li Y., Han B., 2017. Accurate interpolation at receiver positions: A novel method for frequency-domain marine CSEM finite-difference modelling, Pure Appl. Geophys. , 174( 5), 2143– 2160. Google Scholar CrossRef Search ADS   Mackie R.L., Madden T.R., Wannamaker P.E., 1993. Three-dimensional magnetotelluric modeling using difference equations—theory and comparisons to integral equation solutions, Geophysics , 58( 2), 215– 226. Google Scholar CrossRef Search ADS   Martin R., Komatitsch D., Ezziani A., 2009. An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation, Geophysics , 179( 1), 333– 344. Mittet R., 2010. High-order finite-difference simulations of marine CSEM surveys using a correspondence principle for wave and diffusion fields, Geophysics , 75( 1), F33– F50. Google Scholar CrossRef Search ADS   Newman G.A., Alumbaugh D.L., 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences, Geophys. Prospect. , 43( 8), 1021– 1042. Google Scholar CrossRef Search ADS   Oldenburg D.W., Haber E., Shekhtman R., 2013. Three dimensional inversion of multisource time domain electromagnetic data, Geophysics , 78( 1), E47– E57. Google Scholar CrossRef Search ADS   Pan G., Abubakar A., Habashy T.M., 2012. An effective perfectly matched layer design for acoustic fourth-order frequency-domain finite-difference scheme, Geophys. J. Int. , 188, 211– 222. Google Scholar CrossRef Search ADS   Roden J.A., Gedney S.D., 2000. Convolutional PML (CPML): an efficient FDTD implementation of the CFS-PML for arbitrary media, Microw. Opt. Technol. Lett. , 27( 5), 334– 338. Google Scholar CrossRef Search ADS   Sasaki Y., Meju M.A., 2009. Useful characteristics of shallow and deep marine CSEM responses inferred from 3D finite-difference modeling, Geophysics , 74( 5), F67– F76. Google Scholar CrossRef Search ADS   Schwarzbach C., B’orner 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   Smith J.T., 1996. Conservative modeling of 3-D electromagnetic fields, Part I: Properties and error analysis, Geophysics , 61( 5), 1308– 1318. 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   Streich R., Becken M., Ritter O., 2013. Robust processing of noisy land-based controlled-source electromagnetic data, Geophysics , 78( 5), E237– E247. Google Scholar CrossRef Search ADS   Weiss C.J., Constable S., 2006. Mapping thin resistors and hydrocarbons with marine EM methods, Part II—Modeling and analysis in 3D, Geophysics , 71( 6), G321– G332. Google Scholar CrossRef Search ADS   Yee K., 1966. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Trans. Antennas Propag. , 14( 3), 302– 307. Google Scholar CrossRef Search ADS   Zhang Y., 2017. Parallel goal-oriented adaptive finite element modeling for 3D electromagnetic exploration, PhD thesis , University of California, San Diego. Zhdanov M.S., 2010. Electromagnetic geophysics: Notes from the past and the road ahead, Geophysics , 75( 5), 75A49– 75A66. Google Scholar CrossRef Search ADS   APPENDIX A: THE SFD DISCRETIZATION We use the staggered gridding discretization for the governing eq. (3). The SFD discretization around cell (i, j, k) is shown in Fig. A1. For the root grid node (i, j, k) located in the cell centres, the electric field components Ex, Ey and Ez are sampled on the cell faces while the magnetic fields Hx, Hy and Hz are on the edges. Figure A1. View largeDownload slide The 3-D staggered-grid for cell centred at (i, j, k), where Δx(i), Δy(j) or Δz(k) is the cell size along x-, y- and z-directions, respectively. The electric field components Ex, Ey and Ez are sampled at the centre of cell surfaces, while the magnetic field components Hx, Hy and Hz are sampled at the midpoint of cell edges. Figure A1. View largeDownload slide The 3-D staggered-grid for cell centred at (i, j, k), where Δx(i), Δy(j) or Δz(k) is the cell size along x-, y- and z-directions, respectively. The electric field components Ex, Ey and Ez are sampled at the centre of cell surfaces, while the magnetic field components Hx, Hy and Hz are sampled at the midpoint of cell edges. From eq. (3), we obtain   \begin{eqnarray} &&{\frac{1}{\gamma ^h_y}\frac{\partial }{\partial y} \left( \frac{1}{\gamma ^e_x}\frac{\partial E_y^S}{\partial x} - \frac{1}{\gamma ^e_y}\frac{\partial E_x^S}{\partial y} \right) + \frac{1}{\gamma ^h_z}\frac{\partial }{\partial z} \left( \frac{1}{\gamma ^e_x}\frac{\partial E_z^S}{\partial x} - \frac{1}{\gamma ^e_z}\frac{\partial E_x^S}{\partial z} \right)} \nonumber \\ &&{\quad-\, \text{i}\,\omega \mu _0 \sigma ^\ast E_x^S = \text{i}\omega \mu _0 (\sigma ^\ast - \sigma ^{P\ast }) E_x^P, } \end{eqnarray} (A1)  \begin{eqnarray} &&{\frac{1}{\gamma ^h_z}\frac{\partial }{\partial z} \left( \frac{1}{\gamma ^e_y}\frac{\partial E_z^S}{\partial y} - \frac{1}{\gamma ^e_z}\frac{\partial E_y^S}{\partial z} \right) + \frac{1}{\gamma ^h_x}\frac{\partial }{\partial x} \left( \frac{1}{\gamma ^e_y}\frac{\partial E_x^S}{\partial y} - \frac{1}{\gamma ^e_x}\frac{\partial E_y^S}{\partial x} \right)} \nonumber \\ &&{\quad- \text{i}\,\omega \mu _0 \sigma ^\ast E_y^S = \text{i}\omega \mu _0 (\sigma ^\ast - \sigma ^{P\ast }) E_y^P,} \end{eqnarray} (A2)  \begin{eqnarray} &&{\frac{1}{\gamma ^h_x}\frac{\partial }{\partial x} \left( \frac{1}{\gamma ^e_z}\frac{\partial E_x^S}{\partial z} - \frac{1}{\gamma ^e_x}\frac{\partial E_z^S}{\partial x} \right) + \frac{1}{\gamma ^h_y}\frac{\partial }{\partial y} \left( \frac{1}{\gamma ^e_z}\frac{\partial E_y^S}{\partial z} - \frac{1}{\gamma ^e_y}\frac{\partial E_z^S}{\partial y} \right)} \nonumber \\ &&{\quad- \text{i}\,\omega \mu _0 \sigma ^\ast E_z^S = \text{i}\omega \mu _0 (\sigma ^\ast - \sigma ^{P\ast }) E_z^P, } \end{eqnarray} (A3)where $$\gamma ^e_\nu$$ and $$\gamma ^h_\nu$$ (ν = x, y, z) are the complex stretched coordinate factors for electric and magnetic field components, respectively. Following Newman & Alumbaugh (1995) and Streich (2009), the staggered grid FD discretization for the x-, y- and z-components of the electric field evaluated at $$(i+\frac{1}{2},j,k)$$, $$(i,j+\frac{1}{2},k)$$, and $$(i,j,k+\frac{1}{2})$$ can be written, respectively, as   \begin{eqnarray} &&{\frac{1}{\gamma ^{yh}_j}\frac{1}{\Delta y_j}\left[\frac{E^{yS}_{i+1,j+\frac{1}{2},k} - E^{yS}_{i,j+\frac{1}{2},k}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}} - \frac{E^{yS}_{i+1,j-\frac{1}{2},k} - E^{y^S}_{i,j-\frac{1}{2},k}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}}\right]} \nonumber \\ &-& \frac{1}{\gamma ^{yh}_j}\frac{1}{\Delta y_j}\left[\frac{E^{xS}_{i+\frac{1}{2},j+1,k} - E^{xS}_{i+\frac{1}{2},j,k}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}} - \frac{E^{xS}_{i+\frac{1}{2},j,k} - E^{x^S}_{i+\frac{1}{2},j-1,k}}{\gamma ^{ye}_{j-\frac{1}{2}}\Delta y_{j-\frac{1}{2}}}\right] \nonumber \\ &+& \frac{1}{\gamma ^{zh}_k}\frac{1}{\Delta z_k}\left[\frac{E^{zS}_{i+1,j,k+\frac{1}{2}} - E^{zS}_{i,j,k+\frac{1}{2}}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}} - \frac{E^{zS}_{i+1,j,k-\frac{1}{2}} - E^{zS}_{i,j,k-\frac{1}{2}}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}}\right] \nonumber \\ &-& \frac{1}{\gamma ^{zh}_k}\frac{1}{\Delta z_k}\left[\frac{E^{xS}_{i+\frac{1}{2},j,k+1} - E^{x^S}_{i+\frac{1}{2},j,k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}} - \frac{E^{xS}_{i+\frac{1}{2},j,k} - E^{xS}_{i+\frac{1}{2},j,k-1}}{\gamma ^{ze}_{k-\frac{1}{2}}\Delta z_{k-\frac{1}{2}}}\right] \nonumber \\ &-& \text{i}\omega \mu _0\sigma ^\ast _{i+\frac{1}{2},j,k} E^{xS}_{i+\frac{1}{2},j,k} = i\omega \mu _0 (\sigma ^\ast _{i+\frac{1}{2},j,k} - \sigma ^{P\ast }_{i+\frac{1}{2},j,k}) E^{xP}_{i+\frac{1}{2},j,k}, \nonumber \\ \end{eqnarray} (A4)  \begin{eqnarray} &&{\frac{1}{\gamma ^{zh}_k}\frac{1}{\Delta z_k}\left[\frac{E^{zS}_{i,j+1,k+\frac{1}{2}} - E^{zS}_{i,j,k+\frac{1}{2}}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}} - \frac{E^{zS}_{i,j+1,k-\frac{1}{2}} - E^{zS}_{i,j,k-\frac{1}{2}}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}}\right]} \nonumber \\ &-& \frac{1}{\gamma ^{zh}_k}\frac{1}{\Delta z_k}\left[\frac{E^{yS}_{i,j+\frac{1}{2},k+1} - E^{yS}_{i,j+\frac{1}{2},k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}} - \frac{E^{yS}_{i,j+\frac{1}{2},k} - E^{yS}_{i,j+\frac{1}{2},k-1}}{\gamma ^{ze}_{k-\frac{1}{2}}\Delta z_{k-\frac{1}{2}}}\right] \nonumber \\ &+& \frac{1}{\gamma ^{xh}_i}\frac{1}{\Delta x_i}\left[\frac{E^{xS}_{i+\frac{1}{2},j+1,k} - E^{xS}_{i+\frac{1}{2},j,k}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}} - \frac{E^{xS}_{i-\frac{1}{2},j+1,k} - E^{xS}_{i-\frac{1}{2},j,k}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}}\right] \nonumber \\ &-& \frac{1}{\gamma ^{xh}_i}\frac{1}{\Delta x_i}\left[\frac{E^{yS}_{i+1,j+\frac{1}{2},k} - E^{yS}_{i,j+\frac{1}{2},k}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}} - \frac{E^{yS}_{i,j+\frac{1}{2},k} - E^{yS}_{i-1,j+\frac{1}{2},k}}{\gamma ^{xe}_{i-\frac{1}{2}}\Delta x_{i-\frac{1}{2}}}\right] \nonumber \\ &-& \text{i}\omega \mu _0\sigma ^\ast _{i,j+\frac{1}{2},k} E^{yS}_{i,j+\frac{1}{2},k} = i\omega \mu _0 (\sigma ^\ast _{i,j+\frac{1}{2},k} - \sigma ^{P\ast }_{i,j+\frac{1}{2},k}) E^{yP}_{i,j+\frac{1}{2},k}, \nonumber \\ \end{eqnarray} (A5)  \begin{eqnarray} &&{\frac{1}{\gamma ^{xh}_i}\frac{1}{\Delta x_i}\left[\frac{E^{xS}_{i+\frac{1}{2},j,k+1} - E^{xS}_{i+\frac{1}{2},j,k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}} - \frac{E^{xS}_{i-\frac{1}{2},j,k+1} - E^{xS}_{i-\frac{1}{2},j,k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}}\right]} \nonumber \\ &-& \frac{1}{\gamma ^{xh}_i}\frac{1}{\Delta x_i}\left[\frac{E^{zS}_{i+1,j,k+\frac{1}{2}} - E^{zS}_{i,j,k+\frac{1}{2}}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}} - \frac{E^{zS}_{i,j,k+\frac{1}{2}} - E^{zS}_{i-1,j,k+\frac{1}{2}}}{\gamma ^{xe}_{i-\frac{1}{2}}\Delta x_{i-\frac{1}{2}}}\right] \nonumber \\ &+& \frac{1}{\gamma ^{yh}_j}\frac{1}{\Delta y_j}\left[\frac{E^{yS}_{i,j+\frac{1}{2},k+1} - E^{yS}_{i,j+\frac{1}{2},k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}} - \frac{E^{yS}_{i,j-\frac{1}{2},k+1} - E^{yS}_{i,j-\frac{1}{2},k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}}\right] \nonumber \\ &-& \frac{1}{\gamma ^{yh}_j}\frac{1}{\Delta y_j}\left[\frac{E^{zS}_{i,j+1,k+\frac{1}{2}} - E^{zS}_{i,j,k+\frac{1}{2}}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}} - \frac{E^{zS}_{i,j,k+\frac{1}{2}} - E^{zS}_{i,j-1,k+\frac{1}{2}}}{\gamma ^{ye}_{j-\frac{1}{2}}\Delta y_{j-\frac{1}{2}}}\right] \nonumber \\ &-& \text{i}\omega \mu _0\sigma ^\ast _{i,j,k+\frac{1}{2}} E^{zS}_{i,j,k+\frac{1}{2}} = i\omega \mu _0 (\sigma ^\ast _{i,j,k+\frac{1}{2}} - \sigma ^{P\ast }_{i,j,k+\frac{1}{2}}) E^{zP}_{i,j,k+\frac{1}{2}}, \nonumber \\ \end{eqnarray} (A6)where $$\Delta x_{i+\frac{1}{2}}$$, $$\Delta y_{j+\frac{1}{2}}$$ and $$\Delta z_{k+\frac{1}{2}}$$ are the midpoint distances between Δxi and Δxi+1, Δyj and Δyj+1, and Δzk and Δzk+1, respectively. We can obtain an asymmetric system matrix (not Hermitian) equal to that resulting from a finite volume approach given by Weiss & Constable (2006) by multiplying the corresponding cell volume, that is, eq. (A4) by $$(\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}) (\gamma ^{yh}_j\Delta y_j) (\gamma ^{zh}_k\Delta z_k)$$, eq. (A5) by $$(\gamma ^{xh}_i\Delta x_i) (\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}) (\gamma ^{zh}_k\Delta z_k)$$ and eq. (A6) by $$(\gamma ^{xh}_i\Delta x_i) (\gamma ^{yh}_j\Delta y_j) (\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}})$$. We assume model conductivity and permittivity as constant within each cell. As the electric field components are sampled at cell faces where the complex conductivity could be discontinuous, we employ a harmonic averaging method introduced by Smith (1996). For example, the x-, y- and z-staggered complex conductivity is   \begin{eqnarray} \sigma ^\ast _{i+\frac{1}{2},j,k} &=& \frac{2\Delta x_{i+\frac{1}{2}}}{ {\Delta x_i} / {\sigma ^\ast _{i,j,k}} + {\Delta x_{i+1}} / {\sigma ^\ast _{i+1,j,k}} }, \end{eqnarray} (A7)  \begin{eqnarray} \sigma ^\ast _{i,j+\frac{1}{2},k} &=& \frac{2\Delta y_{j+\frac{1}{2}}}{ {\Delta y_j} / {\sigma ^\ast _{i,j,k}} + {\Delta y_{j+1}} / {\sigma ^\ast _{i,j+1,k}} }, \end{eqnarray} (A8)  \begin{eqnarray} \sigma ^\ast _{i,j,k+\frac{1}{2}} &=& \frac{2\Delta z_{k+\frac{1}{2}}}{ {\Delta z_k} / {\sigma ^\ast _{i,j,k}} + {\Delta z_{k+1}} / {\sigma ^\ast _{i,j,k+1}} }, \end{eqnarray} (A9)where $$\Delta x_{i+\frac{1}{2}}$$ is the width of the staggered grid cell that extends between the centres of cell i and i + 1. The harmonic averaging ensures that the current density derived from non-averaged and averaged conductivities cross cell boundaries continuous (Smith 1996; Haber et al. 2000; Streich 2009; Li et al. 2016). For example, the x-, y- and z-staggered complex conductivities are expressed as   \begin{eqnarray} \sigma ^\ast _{i,j,k} E^x_{i,j,k} &=& \sigma ^\ast _{i+1,j,k} E^x_{i+1,j,k} \nonumber \\ &=& \sigma ^\ast _{i+\frac{1}{2},j,k} \frac{\Delta x_i E^x_{i,j,k} + \Delta x_{i+1} E^x_{i+1,j,k}}{\Delta x_i + \Delta x_{i+1}}, \end{eqnarray} (A10)  \begin{eqnarray} \sigma ^\ast _{i,j,k} E^y_{i,j,k} &=& \sigma ^\ast _{i,j+1,k} E^y_{i,j+1,k} \nonumber \\ &=& \sigma ^\ast _{i,j+\frac{1}{2},k} \frac{\Delta y_j E^y_{i,j,k} + \Delta y_{j+1} E^y_{i,j+1,k}}{\Delta y_j + \Delta y_{j+1}}, \end{eqnarray} (A11)  \begin{eqnarray} \sigma ^\ast _{i,j,k} E^z_{i,j,k} &=& \sigma ^\ast _{i,j,k+1} E^z_{i,j,k+1} \nonumber \\ &=& \sigma ^\ast _{i,j,k+\frac{1}{2}} \frac{\Delta z_k E^z_{i,j,k} + \Delta z_{k+1} E^z_{i,j,k+1}}{\Delta z_k + \Delta z_{k+1}}. \end{eqnarray} (A12) © The Authors 2017. 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

Application of the perfectly matched layer in 3-D marine controlled-source electromagnetic modelling

Loading next page...
 
/lp/ou_press/application-of-the-perfectly-matched-layer-in-3-d-marine-controlled-EnsfmlOrev
Publisher
Oxford University Press
Copyright
© The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx382
Publisher site
See Article on Publisher Site

Abstract

SUMMARY In this study, the complex frequency-shifted perfectly matched layer (CFS-PML) in stretching Cartesian coordinates is successfully applied to 3-D frequency-domain marine controlled-source electromagnetic (CSEM) field modelling. The Dirichlet boundary, which is usually used within the traditional framework of EM modelling algorithms, assumes that the electric or magnetic field values are zero at the boundaries. This requires the boundaries to be sufficiently far away from the area of interest. To mitigate the boundary artefacts, a large modelling area may be necessary even though cell sizes are allowed to grow toward the boundaries due to the diffusion of the electromagnetic wave propagation. Compared with the conventional Dirichlet boundary, the PML boundary is preferred as the modelling area of interest could be restricted to the target region and only a few absorbing layers surrounding can effectively depress the artificial boundary effect without losing the numerical accuracy. Furthermore, for joint inversion of seismic and marine CSEM data, if we use the PML for CSEM field simulation instead of the conventional Dirichlet, the modelling area for these two different geophysical data collected from the same survey area could be the same, which is convenient for joint inversion grid matching. We apply the CFS-PML boundary to 3-D marine CSEM modelling by using the staggered finite-difference discretization. Numerical test indicates that the modelling algorithm using the CFS-PML also shows good accuracy compared to the Dirichlet. Furthermore, the modelling algorithm using the CFS-PML shows advantages in computational time and memory saving than that using the Dirichlet boundary. For the 3-D example in this study, the memory saving using the PML is nearly 42  per cent and the time saving is around 48 per cent compared to using the Dirichlet. Controlled source electromagnetics (CSEM), Marine electromagnetics, Numerical modelling 1 INTRODUCTION Marine controlled-source electromagnetic (CSEM) method is widely used for investigating the hydrocarbon and gas hydrate resources (Constable 2010). For frequency-domain marine CSEM survey, a towed electric dipole is usually used to generate the EM fields and the recorded data from receivers located at the seafloor can be used to analyse the resistivity distribution of the seafloor. Nowadays 3-D EM modelling studies are of increasing importance for interpreting the EM data acquired in complex geologic settings (Avdeev 2005; Zhdanov 2010). In the 3-D EM modelling problem, the finite-difference/finite-volume, finite-element and integral equation methods are usually used and these methods are well studied in former publications (e.g. Börner 2010; Kelbert et al. 2014). In this study, the staggered finite-difference (SFD) method is used for modelling the 3-D marine CSEM fields, which can produce accurate and reliable numerical solutions with a fast computation rate (Mackie et al. 1993; Smith 1996; Weiss & Constable 2006; Sasaki & Meju 2009; Mittet 2010). For these traditional frameworks of 3-D EM modelling algorithms, a Dirichlet boundary is usually used, which assumes the electric or magnetic field values are zero at the boundaries. This is a truncated boundary and it requires the boundaries to be sufficiently far away from the area of interest. To mitigate the boundary artefacts, a large modelling area may be necessary though cell sizes are allowed to grow toward the boundaries due to the diffusion of the electromagnetic wave propagation. An alternative way is to use the perfectly matched layer (PML), which is the common choice for seismic wave propagating simulation (e.g. Hu et al. 2007; Pan et al. 2012). By using the PML, the modelling area of interest could be restricted to the target region and only a few surrounding absorbing layers can effectively depress the artificial boundary effect without losing the numerical accuracy. It is worth mentioning that, for joint inversion of marine CSEM and seismic data, if we use the PML for CSEM field simulation instead of the conventional Dirichlet, the modelling area for these two different geophysical data collected from the same survey area could be the same, which is convenient for joint inversion grid matching (Hu et al. 2009). Chen et al. (1997) applied Bérenger’s PML to 2-D transient EM modelling problem using the finite-difference time-domain (FDTD) modelling of lossless media (Bérenger 1994). Schwarzbach et al. (2011) successfully utilized Bérenger’s PML in 3-D frequency-domain CSEM diffusing modelling problem with adaptive higher order finite element method. Similarly, Mittet (2010), de la Kethulle de Ryhove & Mittet (2014) used Bérenger’s PML in 3-D marine magnetotelluric modelling by a finite-difference time-domain algorithm in which Maxwell’s equations are solved in a fictitious-wave domain. The original PML proposed by Bérenger (1994) can be optimized and simply implemented in the form of the complex coordinate stretching along Cartesian coordinates (Chew & Weedon 1994; Chew & Jin 1996). To deduce the asymptotic behaviours at low and high frequencies, a complex frequency-shifted PML (CFS-PML) is also developed which greatly improve the absorbing performance of PML (Kuzuoglu & Mittra 1996; Bérenger 2002). Hu et al. (2017) applied the CFS-PML boundary condition for transient electromagnetic modelling using a fictitious wave domain method. In this study, we present a novel 3-D frequency-domain marine CSEM modelling algorithm, in which CFS-PML is successfully applied. This study can be divided into three parts. We first give an introduction of 3-D frequency-domain CSEM modelling scheme using the SFD grids with the CSF-PML. To avoid the source singularities, the secondary-field approach is used and the primary fields excited by the electric dipole source could be calculated quasi-analytically for the 1-D layered background. Then we compare the performance of CFS-PML with the conventional Dirichlet boundary by numerical analysis. Finally we give the conclusion that the CFS-PML modelling scheme proposed shows advantages in computational time and memory saving compared to using the Dirichlet boundary without losing numerical accuracy. 2 FORMULATIONS 2.1 Governing equations The electric/magnetic fields can be split into a primary part and a secondary (or scattered) part to avoid source singularities (Newman & Alumbaugh 1995; Sasaki & Meju 2009; Streich 2009). The primary fields are excited by the arbitrarily oriented dipole sources and can be computed quasi-analytically following Li & Li (2016). The secondary fields are computed by the frequency-domain SFD method (Yee 1966). Assuming that the time factor is e − iωt where i2 = −1, by using the electric dipole source, the Maxwell’s equations in the frequency domain can be written as   \begin{eqnarray} \nabla \times {{\bf E}}^S - {\rm i}\omega \mu _0 {{\bf H}}^S = 0, \quad \nabla \times {{\bf H}}^S - \sigma ^\ast {{\bf E}}^S = (\sigma ^\ast - \sigma ^{P\ast }) {{\bf E}}^P, \nonumber\\ \end{eqnarray} (1)where superscripts P and S are for the primary and secondary fields, respectively. E and H are the electric and magnetic fields, ω is the angular frequency, μ0 is the magnetic permeability in free space, and the complex conductivity σ* = σ − iωε0 consists of the electric conductivity σ and permittivity ε0 in free space. The secondary-source term is given by (σ* − σP*)EP where σ* = σP* + σS*. After eliminating HS, eq. (1) becomes   \begin{eqnarray} \nabla \times \nabla \times {{\bf E}}^S - {\rm i}\omega \mu _0 \sigma ^\ast {{\bf E}}^S = {\rm i}\omega \mu _0 (\sigma ^\ast - \sigma ^{P\ast }) {{\bf E}}^P. \end{eqnarray} (2)In this study, we use eq. (2) for simulating marine CSEM fields with the frequency range of 0.1–10 Hz. 2.2 Implementation of CFS-PML In the complex stretched coordinate approach (Chew & Weedon 1994; Kuzuoglu & Mittra 1996), the governing equation using CFS-PML can be rewritten as   \begin{eqnarray} \nabla _h \times \nabla _e \times {{\bf E}}^S - {\rm i}\omega \mu _0 \sigma ^\ast {{\bf E}}^S = {\rm i}\omega \mu _0 (\sigma ^\ast - \sigma ^{P\ast }) {{\bf E}}^P. \end{eqnarray} (3)where $$\nabla _e = (\frac{1}{\gamma ^e_x} \partial _x,\frac{1}{\gamma ^e_y} \partial _y,\frac{1}{\gamma ^e_z} \partial _z)$$, $$\nabla _h = (\frac{1}{\gamma ^h_x} \partial _x,\frac{1}{\gamma ^h_y} \partial _y,\frac{1}{\gamma ^h_z} \partial _z)$$, $$\gamma ^e_\nu$$ and $$\gamma ^h_\nu$$ (ν = x, y, z) are the complex stretched coordinate factors (decay factors) for electric and magnetic field components, respectively. SFD discretization of eq. (3) is similar to that in Newman & Alumbaugh (1995) and Streich et al. (2013) (Appendix A). The entire computational domain is divided into Ncell = Nx × Ny × Nz cells, where Nx, Ny, and Nz are the number of cells in the x-, y-, and z- directions, respectively. After discretization, eq. (3) can be written in the matrix form as   \begin{equation} {{\bf K}} {{\bf U}} = {{\bf P}}, \end{equation} (4)where K is a symmetric complex matrix of dimension (3Nx × Ny × Nz)2; the unknown vector U of length 3Nx × Ny × Nz contains the electric field values $$E_x^S$$, $$E_y^S$$ and $$E_z^S$$ for all nodes; and P of length 3Nx × Ny × Nz is the secondary source vector of the right-hand side of eq. (4) given by eq. (3). The entries in K depend on the grid spacing and the frequency-dependent medium properties, and K is up to 13 nonzero entries per line. The 3-D array electric field values $$E_{i+\frac{1}{2},j,k}^{xS}$$, $$E_{i,j+\frac{1}{2},k}^{yS}$$, and $$E_{i,j,k+\frac{1}{2}}^{zS}$$, where the subscripts i = 1, …, Nx, j = 1, …, Ny, and k = 1, …, Nz, can be mapped into a 1-D column array of indices (i + 1/2, j, k) → [(i + 1/2 − 1) × Ny + j] × Nz + k, (i, j + 1/2, k) → [(i − 1) × Ny + j + 1/2] × Nz + k, and (i, j, k + 1/2) → [(i − 1) × Ny + j] × Nz + k + 1/2, respectively. Similarly, we can write the source term in a 1-D column array. The stiffness matrix K is symmetric and highly sparse. Fig. 1 shows an example of the symmetric structure of the stiffness matrix K where a 5 × 5 × 5 staggered grid is used. Figure 1. View largeDownload slide Illustration of the symmetric structure of the stiffness matrix for a 5 × 5 × 5 staggered grid. Note that only the locations of nonzero entries before applying the Dirichlet or PML boundary condition are shown (blue points). nz = 1944 is the number of nonzero entries of the stiffness matrix. Figure 1. View largeDownload slide Illustration of the symmetric structure of the stiffness matrix for a 5 × 5 × 5 staggered grid. Note that only the locations of nonzero entries before applying the Dirichlet or PML boundary condition are shown (blue points). nz = 1944 is the number of nonzero entries of the stiffness matrix. The formed linear equations in (4) given by the discretization of eq. (3) are solved by a multifrontal direct solver MUMPS 5.0.2 parallelized by OpenMP (Amestoy et al. 2001, 2012), which could avoid uncertainties in pre-conditioning and convergence for iterative solutions, especially for low frequencies (Farquharson & Miensopust 2011; Oldenburg et al. 2013). In this section, we will focus on the implementation of CFS-PML in detail. Following Chew & Weedon (1994) and Bérenger (2007b), the PML decay factors are given as   \begin{eqnarray} \gamma ^e_\nu &=& 1 - {\rm i} \eta ^e_\nu / \omega , \quad \nu =x,y,z \end{eqnarray} (5)  \begin{eqnarray} \gamma ^h_\nu &=& 1 - {\rm i} \eta ^h_\nu / \omega , \quad \nu =x,y,z. \end{eqnarray} (6)Absorbing boundaries at the edges of the simulation region may be created by choosing appropriate values of ην. Usually, the decay factor ην varies gradually from 0 at the PML/non-PML interface to its maximum at the outer boundaries of the computational domain to minimize the numerical reflections caused by spatial discretization (Hu et al. 2007). In this study, the PML decay factor ην is determined empirically using the polynomial form of the PML (Hu et al. 2007). For example, the artificial attenuation along the x-direction is defined as (Bérenger 1996; Pan et al. 2012):   \begin{eqnarray} \eta _x = \left\lbrace \begin{array}{rcl}\eta _{\rm {max}} \left( -\frac{x}{L_{\rm {pml}}} \right)^m, & &\quad {{\rm if}\,\quad -L_{\rm {pml}} \le x < 0}\\ {0,} & &\quad {{\rm if}\,\quad 0 \le x \le 1}\\ {\eta _{\rm {max}} \left( \frac{x-1}{L_{\rm {pml}}} \right)^m,} & &\quad {{\rm if}\,\quad 1 < x \le 1+L_{\rm {pml}}} \end{array} \right. \end{eqnarray} (7)where m is either 2 or 3, and this value depends on the size of the computational domain. As pointed out by Bérenger (2002), m = 3 has a better absorption effect on the evanescent region than m = 2, whereas m = 2 has a better absorption effect on the traveling region. We use m = 2 in this study. The length of the interior domain is scaled to 1, Lpml is the thickness of the PML layer on each side, ηmax is the maximum PML decay factor. We use the same definitions for the artificial attenuation parameters in the y- and z-directions. Fig. 2 shows the different PML layers surrounding the original computational area. Note for the computation of the solution inside the PML layers in the x-direction, ηx is positive and ηy, ηz = 0. For PML layers in the y-direction, ηy is positive and ηx, ηz = 0. For PML layers in the z-direction, ηz is positive and both ηx, ηy = 0 (see Table 1). For PML layers in the corners, where we need damping in all directions, ηx, ηy and ηz are all positive (Collino & Tsogka 2001). Figure 2. View largeDownload slide The schematic map showing different PML layers surrounding the 3-D original computational cube. The left plot shows the original computational volume, which is a cuboid without PML layers. The right plot shows the PML layers surrounding the surrounding the original computational volume. The PML layers can be divided into seven different types, which have different values of ηx, ηy and ηz (see Table 1). Figure 2. View largeDownload slide The schematic map showing different PML layers surrounding the 3-D original computational cube. The left plot shows the original computational volume, which is a cuboid without PML layers. The right plot shows the PML layers surrounding the surrounding the original computational volume. The PML layers can be divided into seven different types, which have different values of ηx, ηy and ηz (see Table 1). Table 1. The damping factors ην(ν = x, y, z) for different PML layers. PML layers  PML-x  PML-y  PML-z  PML-x, y  PML-y, z  PML-z, x  PML-x, y, z  ηx  >0  =0  =0  >0  =0  >0  >0  ηy  =0  >0  =0  >0  >0  =0  >0  ηz  =0  =0  >0  =0  >0  >0  >0  PML layers  PML-x  PML-y  PML-z  PML-x, y  PML-y, z  PML-z, x  PML-x, y, z  ηx  >0  =0  =0  >0  =0  >0  >0  ηy  =0  >0  =0  >0  >0  =0  >0  ηz  =0  =0  >0  =0  >0  >0  >0  View Large The complex frequency stretched (CFS) PML, as introduced in Kuzuoglu & Mittra (1996), Roden & Gedney (2000) and Pan et al. (2012), is one of a number of approaches that are suitable to mimic the reflection-free boundaries at the outermost surface of the interior domain. The complex frequency stretched coordinate factors are given by   \begin{eqnarray} \gamma ^e_\nu = \kappa ^e_\nu + \frac{\eta ^e_\nu }{\alpha ^e_\nu + \rm {i}\omega }, \end{eqnarray} (8)  \begin{eqnarray} \gamma ^h_\nu = \kappa ^h_\nu + \frac{\eta ^h_\nu }{\alpha ^h_\nu + \rm {i}\omega }, \end{eqnarray} (9)where ην and κν are positive real and κν is ≥1 (Bérenger 2007a). ην are the PML decay factors, αν are point-wise coefficients related to the CFS-PML (Roden & Gedney 2000). The essence of the real parameter αν is to absorb the evanescent waves and to improve the absorption performance at grazing angles (Bérenger 2002; Komatitsch & Martin 2007; Martin et al. 2009). In this study, the parameter κν is set to be 1, since the effect of these parameters for the absorption is not significant (Bérenger 2002). In the PML region, the conductivity value is equal to that of the outmost cell of the interior domain. This strategy is employed to minimize the numerical reflection error at the interface between the PML domain and the interior domain. We use similar setting for αν as Roden & Gedney (2000). We choose to make αν vary in a linear fashion in their respective PML layer between a maximum value αmax at the beginning (i.e. the entrance) of the PML and zero at its top. The maximum value of αν is set to be 95 per cent of the transmitting frequency. With these settings, the discretized eq. (3) has the same formulas everywhere in space. The only difference between the PML domain and the interior domain equation is the value of the artificial attenuation ην. To reduce the discretization error, the parameter ην are scaled such that they are 0 and 1 at the PML/interior volume interface, respectively, and are maximum at the exterior boundary. We use an optimal maximum attenuation (ηmax) formulation similar to that in Komatitsch & Martin (2007), Martin et al. (2009) and Pan et al. (2012)   \begin{equation} \eta _{\rm {max}} (x,y,z) = c \frac{(m+1)\sigma (x,y,z)}{L_{\rm {pml}}}, \end{equation} (10)where σ(x, y, z) is the conductivity of the cell adjacent to the PML grid, Lpml is the thickness of the PML domain and c a constant that determines the reflection from the PML interface. In this study, the empirical parameter c is chosen to be several times of the skin depth for the background model. 3 NUMERICAL ANALYSIS In this section, we examine the performance of the CFS-PML proposed. Two examples are tested, that are, the 1-D case compared with the quasi-analytical solutions given by Li & Li (2016) and the 3-D case compared with the solutions of adaptive finite-element method (adaptive FEM) given by MARE3DEM code (Zhang 2017). We have written a Fortran 90 code to implement the 3-D SFD modelling scheme for marine CSEM data. The numerical accuracy, time consumption and memory saving are compared between the CFS-PML and the conventional Dirichlet boundary (Li et al. 2017). The numerical test is performed under the Dell Precision Tower 3620 (3.50 GHz CPU Intel@Xeon E3-1240 v5 family including two processors with up to four cores per processor, memory up to 16 G), which is suitable for OpenMP parallel programming. Note the MUMPS 5.0.2 used for factorizing the complex sparse matrix formed by SFD discretization is parallelized by OpenMP in an ‘out-of-core’ environment (Amestoy et al. 2001, 2012). When the ‘out-of-core’ phase is activated and the complete matrix of factors is written to disk and will be read each time a solution phase is requested, therefore the memory requirement can be significantly reduced while not increasing much of the factorization time on a reasonably small number of processors. 3.1 The 1-D case For simplicity, we first use a 1-D canonical reservoir model given by Constable & Weiss (2006), so that our SFD numerical solutions can be easily validated by the results obtained from the quasi-analytical 1-D algorithm (Li & Li 2016). The 1-D canonical reservoir model consists of a sea-water layer with 0.3 Ωm resistivity and 1 km thickness, a sediment layer with 1 Ωm resistivity and a 1000 m thickness (see Fig. 3). The canonical reservoir is buried at a depth of 2 km below the sea surface with 100 Ωm resistivity and a 100 m thickness. Inline (transmitter pointing along y-direction) transmissions are excited at a frequency of 0.25 Hz. The transmitter location is x = 0, y = −5000, z = 950 m, and 51 receivers are located at the seafloor with 200-m intervals from y = −5000 to 5000 m along the towed line. Here, we only consider the electric and magnetic fields in the transmitter-receiver plane at x = 0. Figure 3. View largeDownload slide A 1-D canonical reservoir model similar to Constable & Weiss (2006). There are five layers including air layer, sea-water layer and three seafloor layers. The • and ▾ indicate the transmitter and receivers, respectively. Figure 3. View largeDownload slide A 1-D canonical reservoir model similar to Constable & Weiss (2006). There are five layers including air layer, sea-water layer and three seafloor layers. The • and ▾ indicate the transmitter and receivers, respectively. The computational volume is {(x, y, z): −50 km ≤ x, y, z ≤ 50 km} and it is divided into 76 × 76 × 72 staggered-grids when using the Dirichlet boundary condition (Figs 4a). The x- and y- grid spacings are the same and becomes larger towards the boundaries, and the z-grid spacing becomes larger with depth. The minimum horizontal spacing is set to be 200 m, while the minimum vertical grid spacing is 50 m. For the CFS-PML, the computational volume is {(x, y, z): −6.6 km ≤ x, y ≤ 6.6 km, −8 km ≤ z ≤ 3 km} and a 66 × 66 × 64 staggered-grid including 8 PML layers surrounding the computational area is used (Fig. 4b). Figure 4. View largeDownload slide A 2-D cutaway of the mesh in the transmitter–receiver plane at x = 0 used for forward calculation of the 1-D canonical reservoir model in Fig. 3 using a 3-D SFD code with a PML and Dirichlet boundary condition. (a) The mesh using the conventional Dirichlet boundary in which local mesh for reservoir area bounded by yellow box is enlarged. (b) Local mesh using the CFS-PML boundary. Figure 4. View largeDownload slide A 2-D cutaway of the mesh in the transmitter–receiver plane at x = 0 used for forward calculation of the 1-D canonical reservoir model in Fig. 3 using a 3-D SFD code with a PML and Dirichlet boundary condition. (a) The mesh using the conventional Dirichlet boundary in which local mesh for reservoir area bounded by yellow box is enlarged. (b) Local mesh using the CFS-PML boundary. Fig. 5 shows the 3-D SFD solutions compared to 1-D quasi-analytic solutions for the frequency 0.25 Hz. From Fig. 5 we can see that, the numerical accuracy is acceptable for both the gridding using Dirichlet and CFS-PML boundary. The absolute errors of the amplitude for Ey, Ez and Hx are no more than 2 per cent, 4 per cent and 2 per cent, respectively. The relative errors of the phase for Ey, Ez and Hx are all less than 1°. When using the CFS-PML boundary, less grids are used so that both the memory and computation time saving are significant (Table 2). For this 1-D test, the memory saving is nearly 43  per cent and the time saving is around 50 per cent. Figure 5. View largeDownload slide 3-D numerical solutions against 1-D quasi-analytic solutions for the frequency 0.25 Hz. The y-axis in km is also the source receiver offset. The SFD solutions with the CFS-PML boundary are compared to those with the Dirichlet boundary condition given by Li et al. (2017). (a) Amplitude of the electric field E and electric field H in which the units are VA − 1 m − 2 and m − 2, respectively. (b) Phase in degrees for the electric field E and electric field H. (c) Absolute error of amplitude for the electric field E and electric field H compared to the 1-D analytic solution given by Li & Li (2016). (d) Relative error of phase for the electric field E and electric field H compared to the 1-D analytic solution. The solid line indicates the quasi-analytic results for the 1-D model in Fig. 3, × and ◯ indicate the SFD results using the Dirichlet boundary and the PML boundary, respectively. The colours black, blue and red are for Ey, Ez and Hx, respectively. Figure 5. View largeDownload slide 3-D numerical solutions against 1-D quasi-analytic solutions for the frequency 0.25 Hz. The y-axis in km is also the source receiver offset. The SFD solutions with the CFS-PML boundary are compared to those with the Dirichlet boundary condition given by Li et al. (2017). (a) Amplitude of the electric field E and electric field H in which the units are VA − 1 m − 2 and m − 2, respectively. (b) Phase in degrees for the electric field E and electric field H. (c) Absolute error of amplitude for the electric field E and electric field H compared to the 1-D analytic solution given by Li & Li (2016). (d) Relative error of phase for the electric field E and electric field H compared to the 1-D analytic solution. The solid line indicates the quasi-analytic results for the 1-D model in Fig. 3, × and ◯ indicate the SFD results using the Dirichlet boundary and the PML boundary, respectively. The colours black, blue and red are for Ey, Ez and Hx, respectively. Table 2. Comparison between the CFS-PML boundary and the conventional Dirichlet boundary for the 1-D example. Boundary  Gridding  PML layers  Number of unknowns  Memory used (G)  Total time used (min)  Dirichlet  76 × 76 × 72  –  1247616  ≈7.76  ≈14.15  CFS-PML  66 × 66 × 64  8  836352  ≈4.41  ≈7.04  Boundary  Gridding  PML layers  Number of unknowns  Memory used (G)  Total time used (min)  Dirichlet  76 × 76 × 72  –  1247616  ≈7.76  ≈14.15  CFS-PML  66 × 66 × 64  8  836352  ≈4.41  ≈7.04  View Large Similar to Roden & Gedney (2000) and Hu et al. (2017), the reflection errors at the receiver positions are calculated using   \begin{eqnarray} {\rm Reflection\,\,error} = 20\log _{10} \left| \frac{{\rm abs}({{\bf F}}_i) - {\rm abs}({{\bf F}}^0_i)}{{\rm abs}({{\bf F}}^0_i)} \right|, i=1,\ldots ,N_{R_x} \nonumber\\ \end{eqnarray} (11)where $$N_{R_x}$$ is the number of receivers, F represents the electric E or magnetic field H observed at the receiver locations using the CFS-PML boundary, F0 represents the reference field with boundaries distant enough to avoid any boundary reflections that might interfere with the observed data. In this example, the 1-D quasi-analytical results are taken as the reference fields. Fig. 6 shows the reflection error for both the electric field (Fig. 6a) and the magnetic field (Fig. 6b) at the receiver positions for the 1-D model shown in Fig. 3. Actually, the reflection errors at the receiver positions in eq. (11) indicate the total electric or magnetic field reflection level. The reflection error of both E and H using CFS-PML boundary is low (less than −30 dB noise which corresponds to the relative error of 3 per cent), which is similar to that using the Dirichlet boundary. This indicates that the CFS-PML boundary could effectively depress the artificial boundary effect. Figure 6. View largeDownload slide Reflection error of the PML and Dirichlet boundaries for the 1-D model shown in Fig. 3. The symbols × and ◯ are for the Dirichlet boundary and the PML boundary, respectively. (a) Reflection error for the electric field E, and (b) for the magnetic field H. Figure 6. View largeDownload slide Reflection error of the PML and Dirichlet boundaries for the 1-D model shown in Fig. 3. The symbols × and ◯ are for the Dirichlet boundary and the PML boundary, respectively. (a) Reflection error for the electric field E, and (b) for the magnetic field H. 3.2 The 3-D case In this part, we present a 3-D example similar to that given by Weiss & Constable (2006) (Fig. 7). The 3-D geo-electric model consists of a reservoir represented by a 100-m thick horizontal slab and resistivity 100 Ωm. The reservoir body is buried at a depth of 1 km within a half-space of 1 Ωm water-saturated sediment. The thin slab is a cuboid with 4 km length, 4 km width and 100 m height. The centre of the top surface of the slab corresponds to the origin o of the Cartesian coordinates used. Inline (transmitter pointing along y-direction) transmissions are excited at a frequency of 1 Hz. The transmitter location is x, y = 0, z = 900 m, and 21 receivers are located at the seafloor with 200 m intervals from y = 0 to 5000 m along the towed line. Only the electric and magnetic fields in the transmitter-receiver plane at x = 0 are considered. Figure 7. View largeDownload slide A 3-D canonical reservoir model similar to Weiss & Constable (2006). There are air layer, sea-water layer and a canonical slab of the volume {(x, y, z): −2 km ≤ x, y ≤ 2 km, 2 km ≤ z ≤ 2.1 km} which is embedded into a seafloor half-space. The • and ▾ indicate the transmitter and receivers, respectively. Figure 7. View largeDownload slide A 3-D canonical reservoir model similar to Weiss & Constable (2006). There are air layer, sea-water layer and a canonical slab of the volume {(x, y, z): −2 km ≤ x, y ≤ 2 km, 2 km ≤ z ≤ 2.1 km} which is embedded into a seafloor half-space. The • and ▾ indicate the transmitter and receivers, respectively. The computational volume is {(x, y, z): −50 km ≤ x, y, z ≤ 50 km} is also divided into 78 × 78 × 74 staggered-grids when using the Dirichlet boundary condition (Figs 8a). For the CFS-PML, the computational volume is {(x, y, z): −6.6 ≤ x, y ≤ 6.6 km, −8 km ≤ z ≤ 3 km} and a 68 × 68 × 64 staggered-grid including 8 PML layers surrounding the computational area is used (Fig. 8b). For comparison, the adaptive FEM solutions computed by MARE3DEM are used (Zhang 2017). The computational volume is {(x, y, z): −50 ≤ x, y, z ≤ 50 km} and the total triangular elements are 755 804 after 10 refinement iterations. Figure 8. View largeDownload slide A 2-D cutaway of the mesh in the transmitter-receiver plane at x = 0 used for forward calculation of the 3-D canonical reservoir model in Fig. 7 using a 3-D SFD code with a PML and Dirichlet boundary condition. (a) The mesh using the conventional Dirichlet boundary, in which local mesh for reservoir area bounded by yellow box is enlarged. (b) Local mesh using the CFS-PML boundary. Figure 8. View largeDownload slide A 2-D cutaway of the mesh in the transmitter-receiver plane at x = 0 used for forward calculation of the 3-D canonical reservoir model in Fig. 7 using a 3-D SFD code with a PML and Dirichlet boundary condition. (a) The mesh using the conventional Dirichlet boundary, in which local mesh for reservoir area bounded by yellow box is enlarged. (b) Local mesh using the CFS-PML boundary. Fig. 9 shows the 3-D SFD solutions compared to 3-D adaptive FEM solutions for the frequency 0.25 Hz. Fig. 9 shows a good numerical accuracy for using Dirichlet or CFS-PML boundary. The relative errors of the amplitude for Ey, Ez and Hx are all no more than 1.2 per cent. The absolute errors of the phase are all less than 1°. Both the memory and computation time saving are significant when using the CFS-PML boundary with less grids (Table 3). For the 3-D test, the memory saving using the PML is nearly 42  per cent and the time saving is around 48 per cent compared to using the Dirichlet. In Fig. 9, one can distinguish the 3-D adaptive finite-element results (solid line) and quasi-analytic results for the 1-D background model without the 3-D resistive slab (dashed line). This indicates the effect of the 3-D thin slab is apparent. Figure 9. View largeDownload slide 3-D numerical solutions against 3-D adaptive FEM solutions of MARE3DEM given by Zhang (2017) for the frequency 1 Hz. The SFD solutions with the CFS-PML boundary are compared to those with the Dirichlet boundary condition given by Li et al. (2017). Note that the thin slab is of dimension −4 km to 4 km along y axis. The solid line indicates the adaptive finite-element results, the dash line indicates the quasi-analytic results for the 1-D background model in Fig. 7 without the 3-D resistive slab, × and ◯ indicate the SFD results using the Dirichlet boundary and the PML boundary, respectively. The colours black, blue and red are for Ey, Ez and Hx, respectively. Figure 9. View largeDownload slide 3-D numerical solutions against 3-D adaptive FEM solutions of MARE3DEM given by Zhang (2017) for the frequency 1 Hz. The SFD solutions with the CFS-PML boundary are compared to those with the Dirichlet boundary condition given by Li et al. (2017). Note that the thin slab is of dimension −4 km to 4 km along y axis. The solid line indicates the adaptive finite-element results, the dash line indicates the quasi-analytic results for the 1-D background model in Fig. 7 without the 3-D resistive slab, × and ◯ indicate the SFD results using the Dirichlet boundary and the PML boundary, respectively. The colours black, blue and red are for Ey, Ez and Hx, respectively. Table 3. Comparison between the CFS-PML boundary and the conventional Dirichlet boundary for the 3-D example. Boundary  Gridding  PML layers  Number of unknowns  Memory used (G)  Total time used (min)  Dirichlet  78 × 78 × 74  –  1350648  ≈8.75  ≈15.46  CFS-PML  68 × 68 × 64  8  887808  ≈5.09  ≈8.08  Boundary  Gridding  PML layers  Number of unknowns  Memory used (G)  Total time used (min)  Dirichlet  78 × 78 × 74  –  1350648  ≈8.75  ≈15.46  CFS-PML  68 × 68 × 64  8  887808  ≈5.09  ≈8.08  View Large We notice that the PML appears to have a little higher error compared to the Dirichlet for both the 1-D and 3-D examples. We think that because the modelling area is large enough and the boundaries are far away enough (the modelling area is set to be {(x, y, z): −50 km ≤ x, y, z ≤ 50 km}), no or little boundary reflections appear for the Dirichlet. For the PML, the modelling area for PML is still much smaller than that of Dirichlet. Although we use 10 PML layers which are proved to be sufficient enough for absorbing boundary reflections, the absorbing rate cannot reach 100 per cent. The reflection errors at the receiver positions for the 3-D example are shown in Fig. 10. The 3-D adaptive FEM results, using the Dirichlet boundary, are taken as the reference fields. In this test, the Dirichlet boundary for finite element simulation is far enough to avoid any boundary reflections that might interfere with the observed data. The reflection errors of both E and H using CFS-PML boundary are less than −30 dB and −28 dB, respectively, which are similar to those using the Dirichlet boundary. This indicates that the artificial boundary effect is effectively depressed by the CFS-PML boundary. Figure 10. View largeDownload slide Reflection error of the PML and Dirichlet boundaries for the 3-D model shown in Fig. 7. Note that the thin slab is of dimension −4 km to 4 km along y axis. The symbols × and ◯ are for the Dirichlet boundary and the PML boundary, respectively. (a) Reflection error for the electric field E, and (b) for the magnetic field H. Figure 10. View largeDownload slide Reflection error of the PML and Dirichlet boundaries for the 3-D model shown in Fig. 7. Note that the thin slab is of dimension −4 km to 4 km along y axis. The symbols × and ◯ are for the Dirichlet boundary and the PML boundary, respectively. (a) Reflection error for the electric field E, and (b) for the magnetic field H. 4 CONCLUSIONS In this study, the CFS-PML has been applied to the 3-D marine CSEM forward problem. For the PML boundary, the model area can be restricted to the region of interest and only a few absorbing layers surrounding can effectively depress the artificial boundary effect without losing the numerical accuracy. The next step is to develop an inversion scheme for interpreting real data based on the 3-D marine CSEM modelling algorithm using CFS-PML. Furthermore, for developing a joint 3-D inversion scheme of marine CSEM and seismic data, the proposed 3-D CSEM modelling scheme using CFS-PML could be more convenient than using the conventional Dirichlet boundary condition. For seismic wave simulation, a PML boundary is commonly used. If we use the PML for CSEM field simulation, the modelling area for these two different geophysical data collected from the same survey area could be the same. This could avoid dealing with the grid matching when the CSEM gridding area using Dirichlet boundary is usually larger than the seismic gridding area (Hu et al. 2009). ACKNOWLEDGEMENTS The authors acknowledge the editor Jörg Renner, the associate editor Mark Everett, and two anonymous reviewers for their constructive comments. Wenyi Hu is acknowledged for his fruitful discussions about the implementation of PML. Yuxiang Zhang is thanked for providing the 3-D adaptive finite-element results using MARE3DEM. Gang Li thanks the support from the Sino-German (CSC-DAAD) Postdoc Scholarship Program 2016. This work is funded by the National Natural Science Foundation of China (41774080, 41704075 and 41604063) and the Shandong Provincial Natural Science Foundation, China (ZR2016DQ15 and ZR2016DB31). Finally, Gang Li appreciates the support and help from Marion Jegen and Sebastian Hoelz for his study at GEOMAR. REFERENCES 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   Amestoy P.R., Duff I.S., L’Excellent J.-Y., Robert Y., Rouet F.-H., Uçar B., 2012. On computing inverse entries of a sparse matrix in an out-of-core environment, SIAM J. Sci. Comput. , 34( 4), A1975– A1999. Google Scholar CrossRef Search ADS   Avdeev D.B., 2005. Three-dimensional electromagnetic modelling and inversion from theory to application, Surv. Geophys. , 26( 6), 767– 799. Google Scholar CrossRef Search ADS   Bérenger J.-P., 1994. A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. , 114( 2), 185– 200. Google Scholar CrossRef Search ADS   Bérenger J.-P., 1996. Three-dimensional perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. , 127( 2), 363– 379. Google Scholar CrossRef Search ADS   Bérenger J.-P., 2002. Application of the CFS PML to the absorption of evanescent waves in waveguides, IEEE Microw. Wirel. Compon. Lett. , 12( 6), 218– 220. Google Scholar CrossRef Search ADS   Bérenger J.-P., 2007a. On the Huygens absorbing boundary conditions for electromagnetics, J. Comput. Phys. , 226( 1), 354– 378. Google Scholar CrossRef Search ADS   Bérenger J.-P., 2007b. Perfectly Matched Layer (PML) for Computational Electromagnetics , 1st edn, Morgan & Claypool. Börner R.-U., 2010. Numerical modelling in geo-electromagnetics: advances and challenges, Surv. Geophys. , 31( 2), 225– 245. Google Scholar CrossRef Search ADS   Chen Y.H., Chew W.C., Oristaglio M.L., 1997. Application of perfectly matched layers to the transient modeling of subsurface EM problems, Geophysics , 62( 6), 1730– 1736. Google Scholar CrossRef Search ADS   Chew W.C., Jin J.M., 1996. Perfectly matched layers in the discretized space: an analysis and optimization, Electromagnetics , 16( 4), 325– 340. Google Scholar CrossRef Search ADS   Chew W.C., Weedon W.H., 1994. A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates, Microw. Opt. Technol. Lett. , 7( 13), 599– 604. Google Scholar CrossRef Search ADS   Collino F., Tsogka C., 2001. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media, Geophysics , 66( 1), 294– 307. 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., Weiss C.J., 2006. Mapping thin resistors and hydrocarbons with marine EM methods: insights from 1D modeling, Geophysics , 71( 2), G43– G51. Google Scholar CrossRef Search ADS   de la Kethulle de Ryhove S., Mittet R., 2014. 3D marine magnetotelluric modeling and inversion with the finite-difference time-domain method, Geophysics , 79( 6), E269– E286. Google Scholar CrossRef Search ADS   Farquharson C.G., Miensopust M.P., 2011. Three-dimensional finite-element modelling of magnetotelluric data with a divergence correction, J. Appl. Geophys. , 75( 4), 699– 710. Google Scholar CrossRef Search ADS   Haber E., Ascher U., Aruliah D., Oldenburg D., 2000. Fast simulation of 3D electromagnetic problems using potentials, J. Comput. Phys. , 163( 1), 150– 171. Google Scholar CrossRef Search ADS   Hu W., Abubakar A., Habashy T.M., 2007. Application of the nearly perfectly matched layer in acoustic wave modeling, Geophysics , 72( 5), SM169– SM175. Google Scholar CrossRef Search ADS   Hu W., Abubakar A., Habashy T.M., 2009. Joint electromagnetic and seismic inversion using structural constraints, Geophysics , 74( 6), R99– R109. Google Scholar CrossRef Search ADS   Hu Y., Egbert G., Ji Y., Fang G., 2017. A novel CFS-PML boundary condition for transient electromagnetic simulation using a fictitious wave domain method, Radio Sci. , 52( 1), 118– 131. Google Scholar CrossRef Search ADS   Kelbert A., Kuvshinov A., Velimsky J., Koyama T., Ribaudo J., Sun J., Martinec Z., Weiss C.J., 2014. Global 3-D electromagnetic forward modelling: a benchmark study, Geophys. J. Int. , 197( 2), 785– 814. Google Scholar CrossRef Search ADS   Komatitsch D., Martin R., 2007. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation, Geophysics , 72( 5), SM155– SM167. Google Scholar CrossRef Search ADS   Kuzuoglu M., Mittra R., 1996. Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers, IEEE Microw. Guid. Wave Lett. , 6( 12), 447– 449. Google Scholar CrossRef Search ADS   Li Y., Li G., 2016. Electromagnetic field expressions in the wavenumber domain from both the horizontal and vertical electric dipoles, J. Geophys. Eng. , 13( 4), 505– 515. Google Scholar CrossRef Search ADS   Li G., Zhang L., Han B., 2016. Stable electromagnetic modeling using a multigrid solver on stretching grids: The magnetotelluric example, IEEE Geosci. Remote Sens. Lett. , 13, 334– 338. Li G., Li Y., Han B., 2017. Accurate interpolation at receiver positions: A novel method for frequency-domain marine CSEM finite-difference modelling, Pure Appl. Geophys. , 174( 5), 2143– 2160. Google Scholar CrossRef Search ADS   Mackie R.L., Madden T.R., Wannamaker P.E., 1993. Three-dimensional magnetotelluric modeling using difference equations—theory and comparisons to integral equation solutions, Geophysics , 58( 2), 215– 226. Google Scholar CrossRef Search ADS   Martin R., Komatitsch D., Ezziani A., 2009. An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation, Geophysics , 179( 1), 333– 344. Mittet R., 2010. High-order finite-difference simulations of marine CSEM surveys using a correspondence principle for wave and diffusion fields, Geophysics , 75( 1), F33– F50. Google Scholar CrossRef Search ADS   Newman G.A., Alumbaugh D.L., 1995. Frequency-domain modelling of airborne electromagnetic responses using staggered finite differences, Geophys. Prospect. , 43( 8), 1021– 1042. Google Scholar CrossRef Search ADS   Oldenburg D.W., Haber E., Shekhtman R., 2013. Three dimensional inversion of multisource time domain electromagnetic data, Geophysics , 78( 1), E47– E57. Google Scholar CrossRef Search ADS   Pan G., Abubakar A., Habashy T.M., 2012. An effective perfectly matched layer design for acoustic fourth-order frequency-domain finite-difference scheme, Geophys. J. Int. , 188, 211– 222. Google Scholar CrossRef Search ADS   Roden J.A., Gedney S.D., 2000. Convolutional PML (CPML): an efficient FDTD implementation of the CFS-PML for arbitrary media, Microw. Opt. Technol. Lett. , 27( 5), 334– 338. Google Scholar CrossRef Search ADS   Sasaki Y., Meju M.A., 2009. Useful characteristics of shallow and deep marine CSEM responses inferred from 3D finite-difference modeling, Geophysics , 74( 5), F67– F76. Google Scholar CrossRef Search ADS   Schwarzbach C., B’orner 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   Smith J.T., 1996. Conservative modeling of 3-D electromagnetic fields, Part I: Properties and error analysis, Geophysics , 61( 5), 1308– 1318. 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   Streich R., Becken M., Ritter O., 2013. Robust processing of noisy land-based controlled-source electromagnetic data, Geophysics , 78( 5), E237– E247. Google Scholar CrossRef Search ADS   Weiss C.J., Constable S., 2006. Mapping thin resistors and hydrocarbons with marine EM methods, Part II—Modeling and analysis in 3D, Geophysics , 71( 6), G321– G332. Google Scholar CrossRef Search ADS   Yee K., 1966. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media, IEEE Trans. Antennas Propag. , 14( 3), 302– 307. Google Scholar CrossRef Search ADS   Zhang Y., 2017. Parallel goal-oriented adaptive finite element modeling for 3D electromagnetic exploration, PhD thesis , University of California, San Diego. Zhdanov M.S., 2010. Electromagnetic geophysics: Notes from the past and the road ahead, Geophysics , 75( 5), 75A49– 75A66. Google Scholar CrossRef Search ADS   APPENDIX A: THE SFD DISCRETIZATION We use the staggered gridding discretization for the governing eq. (3). The SFD discretization around cell (i, j, k) is shown in Fig. A1. For the root grid node (i, j, k) located in the cell centres, the electric field components Ex, Ey and Ez are sampled on the cell faces while the magnetic fields Hx, Hy and Hz are on the edges. Figure A1. View largeDownload slide The 3-D staggered-grid for cell centred at (i, j, k), where Δx(i), Δy(j) or Δz(k) is the cell size along x-, y- and z-directions, respectively. The electric field components Ex, Ey and Ez are sampled at the centre of cell surfaces, while the magnetic field components Hx, Hy and Hz are sampled at the midpoint of cell edges. Figure A1. View largeDownload slide The 3-D staggered-grid for cell centred at (i, j, k), where Δx(i), Δy(j) or Δz(k) is the cell size along x-, y- and z-directions, respectively. The electric field components Ex, Ey and Ez are sampled at the centre of cell surfaces, while the magnetic field components Hx, Hy and Hz are sampled at the midpoint of cell edges. From eq. (3), we obtain   \begin{eqnarray} &&{\frac{1}{\gamma ^h_y}\frac{\partial }{\partial y} \left( \frac{1}{\gamma ^e_x}\frac{\partial E_y^S}{\partial x} - \frac{1}{\gamma ^e_y}\frac{\partial E_x^S}{\partial y} \right) + \frac{1}{\gamma ^h_z}\frac{\partial }{\partial z} \left( \frac{1}{\gamma ^e_x}\frac{\partial E_z^S}{\partial x} - \frac{1}{\gamma ^e_z}\frac{\partial E_x^S}{\partial z} \right)} \nonumber \\ &&{\quad-\, \text{i}\,\omega \mu _0 \sigma ^\ast E_x^S = \text{i}\omega \mu _0 (\sigma ^\ast - \sigma ^{P\ast }) E_x^P, } \end{eqnarray} (A1)  \begin{eqnarray} &&{\frac{1}{\gamma ^h_z}\frac{\partial }{\partial z} \left( \frac{1}{\gamma ^e_y}\frac{\partial E_z^S}{\partial y} - \frac{1}{\gamma ^e_z}\frac{\partial E_y^S}{\partial z} \right) + \frac{1}{\gamma ^h_x}\frac{\partial }{\partial x} \left( \frac{1}{\gamma ^e_y}\frac{\partial E_x^S}{\partial y} - \frac{1}{\gamma ^e_x}\frac{\partial E_y^S}{\partial x} \right)} \nonumber \\ &&{\quad- \text{i}\,\omega \mu _0 \sigma ^\ast E_y^S = \text{i}\omega \mu _0 (\sigma ^\ast - \sigma ^{P\ast }) E_y^P,} \end{eqnarray} (A2)  \begin{eqnarray} &&{\frac{1}{\gamma ^h_x}\frac{\partial }{\partial x} \left( \frac{1}{\gamma ^e_z}\frac{\partial E_x^S}{\partial z} - \frac{1}{\gamma ^e_x}\frac{\partial E_z^S}{\partial x} \right) + \frac{1}{\gamma ^h_y}\frac{\partial }{\partial y} \left( \frac{1}{\gamma ^e_z}\frac{\partial E_y^S}{\partial z} - \frac{1}{\gamma ^e_y}\frac{\partial E_z^S}{\partial y} \right)} \nonumber \\ &&{\quad- \text{i}\,\omega \mu _0 \sigma ^\ast E_z^S = \text{i}\omega \mu _0 (\sigma ^\ast - \sigma ^{P\ast }) E_z^P, } \end{eqnarray} (A3)where $$\gamma ^e_\nu$$ and $$\gamma ^h_\nu$$ (ν = x, y, z) are the complex stretched coordinate factors for electric and magnetic field components, respectively. Following Newman & Alumbaugh (1995) and Streich (2009), the staggered grid FD discretization for the x-, y- and z-components of the electric field evaluated at $$(i+\frac{1}{2},j,k)$$, $$(i,j+\frac{1}{2},k)$$, and $$(i,j,k+\frac{1}{2})$$ can be written, respectively, as   \begin{eqnarray} &&{\frac{1}{\gamma ^{yh}_j}\frac{1}{\Delta y_j}\left[\frac{E^{yS}_{i+1,j+\frac{1}{2},k} - E^{yS}_{i,j+\frac{1}{2},k}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}} - \frac{E^{yS}_{i+1,j-\frac{1}{2},k} - E^{y^S}_{i,j-\frac{1}{2},k}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}}\right]} \nonumber \\ &-& \frac{1}{\gamma ^{yh}_j}\frac{1}{\Delta y_j}\left[\frac{E^{xS}_{i+\frac{1}{2},j+1,k} - E^{xS}_{i+\frac{1}{2},j,k}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}} - \frac{E^{xS}_{i+\frac{1}{2},j,k} - E^{x^S}_{i+\frac{1}{2},j-1,k}}{\gamma ^{ye}_{j-\frac{1}{2}}\Delta y_{j-\frac{1}{2}}}\right] \nonumber \\ &+& \frac{1}{\gamma ^{zh}_k}\frac{1}{\Delta z_k}\left[\frac{E^{zS}_{i+1,j,k+\frac{1}{2}} - E^{zS}_{i,j,k+\frac{1}{2}}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}} - \frac{E^{zS}_{i+1,j,k-\frac{1}{2}} - E^{zS}_{i,j,k-\frac{1}{2}}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}}\right] \nonumber \\ &-& \frac{1}{\gamma ^{zh}_k}\frac{1}{\Delta z_k}\left[\frac{E^{xS}_{i+\frac{1}{2},j,k+1} - E^{x^S}_{i+\frac{1}{2},j,k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}} - \frac{E^{xS}_{i+\frac{1}{2},j,k} - E^{xS}_{i+\frac{1}{2},j,k-1}}{\gamma ^{ze}_{k-\frac{1}{2}}\Delta z_{k-\frac{1}{2}}}\right] \nonumber \\ &-& \text{i}\omega \mu _0\sigma ^\ast _{i+\frac{1}{2},j,k} E^{xS}_{i+\frac{1}{2},j,k} = i\omega \mu _0 (\sigma ^\ast _{i+\frac{1}{2},j,k} - \sigma ^{P\ast }_{i+\frac{1}{2},j,k}) E^{xP}_{i+\frac{1}{2},j,k}, \nonumber \\ \end{eqnarray} (A4)  \begin{eqnarray} &&{\frac{1}{\gamma ^{zh}_k}\frac{1}{\Delta z_k}\left[\frac{E^{zS}_{i,j+1,k+\frac{1}{2}} - E^{zS}_{i,j,k+\frac{1}{2}}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}} - \frac{E^{zS}_{i,j+1,k-\frac{1}{2}} - E^{zS}_{i,j,k-\frac{1}{2}}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}}\right]} \nonumber \\ &-& \frac{1}{\gamma ^{zh}_k}\frac{1}{\Delta z_k}\left[\frac{E^{yS}_{i,j+\frac{1}{2},k+1} - E^{yS}_{i,j+\frac{1}{2},k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}} - \frac{E^{yS}_{i,j+\frac{1}{2},k} - E^{yS}_{i,j+\frac{1}{2},k-1}}{\gamma ^{ze}_{k-\frac{1}{2}}\Delta z_{k-\frac{1}{2}}}\right] \nonumber \\ &+& \frac{1}{\gamma ^{xh}_i}\frac{1}{\Delta x_i}\left[\frac{E^{xS}_{i+\frac{1}{2},j+1,k} - E^{xS}_{i+\frac{1}{2},j,k}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}} - \frac{E^{xS}_{i-\frac{1}{2},j+1,k} - E^{xS}_{i-\frac{1}{2},j,k}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}}\right] \nonumber \\ &-& \frac{1}{\gamma ^{xh}_i}\frac{1}{\Delta x_i}\left[\frac{E^{yS}_{i+1,j+\frac{1}{2},k} - E^{yS}_{i,j+\frac{1}{2},k}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}} - \frac{E^{yS}_{i,j+\frac{1}{2},k} - E^{yS}_{i-1,j+\frac{1}{2},k}}{\gamma ^{xe}_{i-\frac{1}{2}}\Delta x_{i-\frac{1}{2}}}\right] \nonumber \\ &-& \text{i}\omega \mu _0\sigma ^\ast _{i,j+\frac{1}{2},k} E^{yS}_{i,j+\frac{1}{2},k} = i\omega \mu _0 (\sigma ^\ast _{i,j+\frac{1}{2},k} - \sigma ^{P\ast }_{i,j+\frac{1}{2},k}) E^{yP}_{i,j+\frac{1}{2},k}, \nonumber \\ \end{eqnarray} (A5)  \begin{eqnarray} &&{\frac{1}{\gamma ^{xh}_i}\frac{1}{\Delta x_i}\left[\frac{E^{xS}_{i+\frac{1}{2},j,k+1} - E^{xS}_{i+\frac{1}{2},j,k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}} - \frac{E^{xS}_{i-\frac{1}{2},j,k+1} - E^{xS}_{i-\frac{1}{2},j,k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}}\right]} \nonumber \\ &-& \frac{1}{\gamma ^{xh}_i}\frac{1}{\Delta x_i}\left[\frac{E^{zS}_{i+1,j,k+\frac{1}{2}} - E^{zS}_{i,j,k+\frac{1}{2}}}{\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}} - \frac{E^{zS}_{i,j,k+\frac{1}{2}} - E^{zS}_{i-1,j,k+\frac{1}{2}}}{\gamma ^{xe}_{i-\frac{1}{2}}\Delta x_{i-\frac{1}{2}}}\right] \nonumber \\ &+& \frac{1}{\gamma ^{yh}_j}\frac{1}{\Delta y_j}\left[\frac{E^{yS}_{i,j+\frac{1}{2},k+1} - E^{yS}_{i,j+\frac{1}{2},k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}} - \frac{E^{yS}_{i,j-\frac{1}{2},k+1} - E^{yS}_{i,j-\frac{1}{2},k}}{\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}}}\right] \nonumber \\ &-& \frac{1}{\gamma ^{yh}_j}\frac{1}{\Delta y_j}\left[\frac{E^{zS}_{i,j+1,k+\frac{1}{2}} - E^{zS}_{i,j,k+\frac{1}{2}}}{\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}} - \frac{E^{zS}_{i,j,k+\frac{1}{2}} - E^{zS}_{i,j-1,k+\frac{1}{2}}}{\gamma ^{ye}_{j-\frac{1}{2}}\Delta y_{j-\frac{1}{2}}}\right] \nonumber \\ &-& \text{i}\omega \mu _0\sigma ^\ast _{i,j,k+\frac{1}{2}} E^{zS}_{i,j,k+\frac{1}{2}} = i\omega \mu _0 (\sigma ^\ast _{i,j,k+\frac{1}{2}} - \sigma ^{P\ast }_{i,j,k+\frac{1}{2}}) E^{zP}_{i,j,k+\frac{1}{2}}, \nonumber \\ \end{eqnarray} (A6)where $$\Delta x_{i+\frac{1}{2}}$$, $$\Delta y_{j+\frac{1}{2}}$$ and $$\Delta z_{k+\frac{1}{2}}$$ are the midpoint distances between Δxi and Δxi+1, Δyj and Δyj+1, and Δzk and Δzk+1, respectively. We can obtain an asymmetric system matrix (not Hermitian) equal to that resulting from a finite volume approach given by Weiss & Constable (2006) by multiplying the corresponding cell volume, that is, eq. (A4) by $$(\gamma ^{xe}_{i+\frac{1}{2}}\Delta x_{i+\frac{1}{2}}) (\gamma ^{yh}_j\Delta y_j) (\gamma ^{zh}_k\Delta z_k)$$, eq. (A5) by $$(\gamma ^{xh}_i\Delta x_i) (\gamma ^{ye}_{j+\frac{1}{2}}\Delta y_{j+\frac{1}{2}}) (\gamma ^{zh}_k\Delta z_k)$$ and eq. (A6) by $$(\gamma ^{xh}_i\Delta x_i) (\gamma ^{yh}_j\Delta y_j) (\gamma ^{ze}_{k+\frac{1}{2}}\Delta z_{k+\frac{1}{2}})$$. We assume model conductivity and permittivity as constant within each cell. As the electric field components are sampled at cell faces where the complex conductivity could be discontinuous, we employ a harmonic averaging method introduced by Smith (1996). For example, the x-, y- and z-staggered complex conductivity is   \begin{eqnarray} \sigma ^\ast _{i+\frac{1}{2},j,k} &=& \frac{2\Delta x_{i+\frac{1}{2}}}{ {\Delta x_i} / {\sigma ^\ast _{i,j,k}} + {\Delta x_{i+1}} / {\sigma ^\ast _{i+1,j,k}} }, \end{eqnarray} (A7)  \begin{eqnarray} \sigma ^\ast _{i,j+\frac{1}{2},k} &=& \frac{2\Delta y_{j+\frac{1}{2}}}{ {\Delta y_j} / {\sigma ^\ast _{i,j,k}} + {\Delta y_{j+1}} / {\sigma ^\ast _{i,j+1,k}} }, \end{eqnarray} (A8)  \begin{eqnarray} \sigma ^\ast _{i,j,k+\frac{1}{2}} &=& \frac{2\Delta z_{k+\frac{1}{2}}}{ {\Delta z_k} / {\sigma ^\ast _{i,j,k}} + {\Delta z_{k+1}} / {\sigma ^\ast _{i,j,k+1}} }, \end{eqnarray} (A9)where $$\Delta x_{i+\frac{1}{2}}$$ is the width of the staggered grid cell that extends between the centres of cell i and i + 1. The harmonic averaging ensures that the current density derived from non-averaged and averaged conductivities cross cell boundaries continuous (Smith 1996; Haber et al. 2000; Streich 2009; Li et al. 2016). For example, the x-, y- and z-staggered complex conductivities are expressed as   \begin{eqnarray} \sigma ^\ast _{i,j,k} E^x_{i,j,k} &=& \sigma ^\ast _{i+1,j,k} E^x_{i+1,j,k} \nonumber \\ &=& \sigma ^\ast _{i+\frac{1}{2},j,k} \frac{\Delta x_i E^x_{i,j,k} + \Delta x_{i+1} E^x_{i+1,j,k}}{\Delta x_i + \Delta x_{i+1}}, \end{eqnarray} (A10)  \begin{eqnarray} \sigma ^\ast _{i,j,k} E^y_{i,j,k} &=& \sigma ^\ast _{i,j+1,k} E^y_{i,j+1,k} \nonumber \\ &=& \sigma ^\ast _{i,j+\frac{1}{2},k} \frac{\Delta y_j E^y_{i,j,k} + \Delta y_{j+1} E^y_{i,j+1,k}}{\Delta y_j + \Delta y_{j+1}}, \end{eqnarray} (A11)  \begin{eqnarray} \sigma ^\ast _{i,j,k} E^z_{i,j,k} &=& \sigma ^\ast _{i,j,k+1} E^z_{i,j,k+1} \nonumber \\ &=& \sigma ^\ast _{i,j,k+\frac{1}{2}} \frac{\Delta z_k E^z_{i,j,k} + \Delta z_{k+1} E^z_{i,j,k+1}}{\Delta z_k + \Delta z_{k+1}}. \end{eqnarray} (A12) © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

There are no references for this article.

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
discover and read the research
that matters to you.

Enjoy affordable access to
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.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off