# Trapezoid-Grid Finite-Difference Time-Domain Method for 3D Seismic Wavefield Modeling Using CPML Absorbing Boundary Condition Bangyu Wu, et al.

Jan 12, 2022

## 1 Introduction

Reverse time migration (RTM) (Baysal et al., 1983; Xuan et al., 2014; Qu et al., 2015; Xu et al., 2021a; Du et al., 2021) and full-waveform inversion (FWI) (Tarantola, 1984; Virieux and Operto, 2009; Cai and Zhang, 2015; Xia et al., 2017; Jia et al., 2019) play a fundamental role in geophysical exploration. Since forward modeling of the wave equation consumes most computational time in the RTM and FWI processes (Jing et al., 2019; Xu et al., 2021b; Liu et al., 2021), how to achieve the improvement of efficiency and the reduction of memory usage without a significant decrease of accuracy for 3D large-scale modeling is a key problem of seismic modeling. The finite-difference (FD) method has been regarded as one of the most popular wave modeling methods for its easy implementation and high-computational efficiency (Antunes et al., 2014; Abreu et al., 2015; Xu and Gao, 2018; Robertsson and Blanch, 2020). However, the numerical dispersion of the traditional FD method leads to the use of fine grids or high-order operators (Dablain, 1986; Liu and Sen, 2011b), which inevitably affects the efficiency of simulation.

The conventional FD method literally adopts a weighted summation of neighboring grid points’ values to estimate the derivative for a designated grid point (Zhou et al., 2021), where the grid size (h) is fixed and the FD coefficients are calculated by Taylor expansion. In this way, the approximation error ϵ can be expressed as (Liu and Sen, 2011a; Wu et al., 2019b) follows:

where 2M is the length of the FD operator, h is the spatial interval, f is the frequency, and v is the seismic wave velocity. Considering that λ = v/f is the wavelength and G = λ/h is the number of grid points per wavelength (NPPW), we can rewrite Eq. 1 as follows:

Eq. 2 indicates that the modeling accuracy of the conventional FD method is proportional to G and M. Because the seismic wave velocity is varying in different positions, the wavelength is short in low-velocity regions and long in high-velocity regions (Liu, 2020). Therefore, a part of computing resources is wasted in the high-velocity regions for the fixed spatial interval and FD order. With respect to this problem, there are two kinds of techniques corresponding to the different understanding of Eq. 2. The first one is the variable-operator FD method (Liu and Sen, 2011a), which adopts the long and short FD stencils in the low- and high-velocity region, respectively. For the scheme in Liu and Sen (2011a), the variable-length FD stencils are designed by approaching the dispersion relation in the time–space domain, and Liu (2020) subsequently optimizes their FD coefficients.

The second one is the variable-grid FD method, which adopts different gird sizes in different regions and can efficiently reduce the oversampling in the high-velocity region. The key problem of the variable-grid FD method is the processing of the transition area between the fine grids and the coarse grids. The variable-grid FD method based on interpolation (Hayashi and Burns, 2005; Pasalic and McGarry, 2010) is the easiest one, in which the lacking information in the transition area is completed by interpolation. However, the resulting artificial reflection in the transition area and the possible instability make it inefficient for high-accuracy seismic wave simulation. Another variable-grid FD method adopts irregular FD coefficients to process the transition region (Huang and Dong, 2009; Liu et al., 2014), which can significantly avoid the artificial reflection and improve the stability. The disadvantage of this type of variable-grid method is the additional computing cost brought by calculating irregular FD coefficients.

The trapezoid-grid FD method (Chen and Xu, 2012; Gao et al., 2018; Wu et al., 2018, Wu et al., 2019a, Wu et al., 2019b) is one of the practical variable-grid methods. It uses the trapezoid-grid mesh to fit the trend of velocity increasing with depth, which can effectively reduce the number of required grid points. Meanwhile, the use of trapezoid coordinate transformation can avoid the difficulty of processing ununiform grids in the physical Cartesian coordinate system. On the other hand, the significant reduction of memory requirement of trapezoid-grid FDTD can improve the easy implementation of GPU calculation (Fujii et al., 2013; Li et al., 2016). The existing research on trapezoid-grid FDTD methods mainly focuses on 2D wavefield modeling. Therefore, it is essential to expand trapezoid-grid FDTD from 2D to 3D for realistic seismic exploration research.

In this article, we propose a 3D trapezoid-grid FDTD method for acoustic wave modeling. First, we design the 3D trapezoid coordinate transformation and derive the 3D acoustic equation in the trapezoid coordinate system. Second, to reduce the artificial boundary reflection (Ma et al., 2018, Ma et al., 2019), we apply the corresponding trapezoid-grid convolutional perfectly matched layer (CPML) absorbing boundary condition. Third, stability analysis is given to generate stable modeling results. We then test our proposed method on the 3D homogenous model and the SEG/EAGE overthrust model and compare the efficiency and accuracy of the trapezoid-grid FDTD method with the uniform-grid FDTD method. Finally, conclusions are shown in the last section.

## 2 Methods

### 2.1 3D Trapezoid Coordinate System

where (x0, y0, z0) is the Cartesian coordinate system, and (x, y, z) is the defined trapezoid coordinate system. In Eq. 3, α and β are central horizontal positions of the 3D trapezoid mesh, and γ is the scaling parameter for lateral coordinates. The velocity-related function g(z) is the sampling function for z0-axis, which should be first- and second-order continuous for deriving 3D wave equations in the trapezoid coordinate system. The discrete points of g(z) are given by the following recursion:

$gz+Δz=gz+vmingzf0N0,(4b)$

where f0 is the dominant frequency of the source term, N0 is the preferred NPPW and is related to the accuracy of the adopted FD scheme, and vmin(g(z)) is the selected minimum velocity at depth g(z) in the physical Cartesian coordinate system and is smoothed by solving a local polynomial fitting problem with the constraint that vmin(g(z)) should not be greater than the model minimum velocity at depth g(z). The central value of each local polynomial corresponds to a value of vmin(g(z)). In particular, we usually set the order of the local polynomial as three. Such vmin(g(z)) can lead to a smooth sampling function g(z) for discontinuous velocity variation while satisfying the required number of points per wavelength in the z0-direction.

If the grid sizes for the trapezoid-grid FDTD in the trapezoid coordinate system are defined as Δx, Δy, and Δz, then the corresponding grid sizes in the physical Cartesian system can be described as

In our work, γ and g(z) are determined adaptively according to the model velocity. By selecting γ such that Δx0(z) and Δy0(z) are always smaller than or equal to Δz0(z), and a variable-grid mesh adaptive to the velocity model can be achieved in the physical Cartesian coordinate system. Figure 1 shows the schematic of the 3D trapezoid coordinate transform. Figure 1A shows the trapezoid-grid mesh in the Cartesian coordinate system, while Figure 1B shows the corresponding uniform grid mesh in the transformed trapezoid coordinate system. In particular, the two gray regions in Figure 1 represent the same physical region.

FIGURE 1. Schematic of the 3D trapezoid coordinate transformation: (A) the trapezoid-grid mesh in the Cartesian coordinate system; (B) the uniform grid mesh in the transformed trapezoid coordinate system. The two gray regions in A and B represent the same physical region.

### 2.2 3D Acoustic Equation with CPML Absorbing Boundary Condition in the Trapezoid Coordinate System

According to the theory of Pasalic and McGarry (2010), the time-domain-discretization form of the 3D isotropic acoustic equation with the CPML absorbing boundary condition in the Cartesian coordinate system can be described as

$1v2uj+1−2uj+uj−1Δt2−∂2uj∂x02−∂2uj∂y02−∂2uj∂z02−∂ψx0j∂x0−∂ψy0j∂y0−∂ψz0j∂z0−ζx0j−ζy0j−ζz0j=ftjδx0−x0sδy0−y0sδz0−z0s;(6a)$
$ψτ0j+1=aτ0ψτ0j+bτ0∂uj+1∂τ0;(6b)$
$ζτ0j+1=aτ0ζτ0j+bτ0∂2uj+1∂τ02+∂ψτ0j+1∂τ0;(6c)$
$aτ0=e−στ0+ατ0Δt;(6d)$
$bτ0=στ0στ0+ατ0aτ0−1,τ0∈x0,y0,z0,(6e)$

where uj = u(x0, y0, z0, tj) represents the scalar wavefield at the jth time step in the Cartesian coordinate system; v is the velocity; Δt is the time interval; f(t) is the source term;

$(x0s,y0s,z0s)$

is the position of source;

$στ0=3vmax2Lτ0τ̄0Lτ02$

ln

$1R$

, where R denotes the designated theoretical boundary reflection coefficient, vmax is the maximum velocity of the model,

$Lτ0$

is the thickness of CPML absorbing boundary along the τ0 direction, and

$τ̄0$

denotes the distance to the inner area in the τ0 direction;

$ατ0=αmax1−τ̄0Lτ0$

and αmax = πf0.

In order to derive the acoustic equation in the trapezoid coordinate system, we first need to transform the derivatives in the Cartesian coordinate system into the derivatives in the trapezoid coordinate system. Based on the definition of the trapezoid coordinate system in Eq. 3 and the derivation rule of the composite function, the relationships of first- and second-order derivatives in the two coordinate systems can be given as

$∂∂x0=11+γgz∂∂x;(7a)$
$∂∂y0=11+γgz∂∂y;(7b)$
$∂∂z0=−γx1+γgz∂∂x−γy1+γgz∂∂y+1g′z∂∂z;(7c)$
$∂2∂x02=11+γgz2∂2∂x2;(7d)$
$∂2∂y02=11+γgz2∂2∂y2;(7e)$
$∂2∂z02=2γ2x1+γgz2∂∂x+γ2x21+γgz2∂2∂x2+2γ2xy1+γgz2∂2∂x∂y−2γx1+γgzg′z∂2∂x∂z+2γ2y1+γgz2∂∂y+γ2y21+γgz2∂2∂y2−2γy1+γgzg′z∂2∂y∂z−g″zg′z3∂∂z+1g′z2∂2∂z2.(7f)$

Since the mixed spatial derivatives in Eq. 7f are hard to discrete directly with the FD method, to transform the mixed spatial derivatives into non-mixed spatial derivatives, we define three rotation transformations in the trapezoid coordinate system as

$xz=cosθ1−sinθ1sinθ1cosθ1x̂ẑ;(8a)$
$yz=cosθ2−sinθ2sinθ2cosθ2ỹz̃;(8b)$
$xy=cosθ3−sinθ3sinθ3cosθ3x̄ȳ,(8c)$

where

$x̂$

and

$ẑ$

are axes along diagonals in the (x, z) planes,

$ỹ$

and

$z̃$

are axes along diagonals in the (y, z) planes, and

$x̄$

and

$ȳ$

are axes along diagonals in the (x, y) planes. θ1 is the angle between x and

$x̂$

axes, θ2 is the angle between y and

$ỹ$

axes, and θ3 is the angle between x and

$x̄$

axes. A schematic of the coordinate transform in Eq. 8a is shown in Figure 2.

FIGURE 2. Schematic of the rotation transformation in the trapezoid coordinate system for transforming mixed spatial derivatives into non-mixed spatial derivatives.

By using Eqs. 8a–c, the mixed spatial derivatives in Eq. 7f can be transformed as

$∂2∂x∂z=12sin2θ1∂2∂x̂2−∂2∂ẑ2;(9a)$
$∂2∂y∂z=12sin2θ2∂2∂ỹ2−∂2∂z̃2;(9b)$
$∂2∂x∂y=12sin2θ3∂2∂x̄2−∂2∂ȳ2.(9c)$

For simplicity, we usually use equal grid sizes in the trapezoid coordinate system (Δx = Δy = Δz = Δ), which means

$θ1=θ2=θ3=π4$

. By substituting Eq. 7 and Eq. 9 into Eq. 6, we get the time-domain-discretization form of the 3D acoustic equation with the CPML absorbing boundary condition in the trapezoid coordinate system as

$1v2uj+1−2uj+uj−1Δt2−γ2x2+11+γgz2∂2uj∂x2−γ2y2+11+γgz2∂2uj∂y2−2γ2x1+γgz2∂uj∂x−2γ2y1+γgz2∂uj∂y+γx1+γgzg′z∂2uj∂x̂2−∂2uj∂ẑ2+γy1+γgzg′z∂2uj∂ỹ2−∂2uj∂z̃2−γ2xy1+γgz2∂2uj∂x̄2−∂2uj∂ȳ2+g″zg′z3∂uj∂z−1g′z2∂2uj∂z2−11+γgz∂ψxj∂x−11+γgz∂ψyj∂y+γx1+γgz∂ψzj∂z+γy1+γgz∂ψzj∂y−1g′z∂ψzj∂z−ζxj−ζyj−ζzj=ftjx−xsy−ysz−zs;(10a)$
$ψxj+1=axψxj+bx11+γgz∂uj+1∂x;(10b)$
$ψyj+1=ayψyj+by11+γgz∂uj+1∂y;(10c)$
$ψzj+1=azψzj+bz−γx1+γgz∂uj+1∂x−γy1+γgz∂uj+1∂y+1g′z∂uj+1∂z;(10d)$
$ζxj+1=axζxj+bx11+γgz2∂2uj+1∂x2+11+γgz∂ψxj+1∂x;(10e)$
$ζyj+1=ayζyj+by11+γgz2∂2uj+1∂y2+11+γgz∂ψyj+1∂y;(10f)$
$ζzj+1=azζzj+bzγ2x21+γgz2∂2uj+1∂x2+γ2y21+γgz2∂2uj+1∂y2+2γ2x1+γgz2∂uj+1∂x+2γ2y1+γgz2∂uj+1∂y−γx1+γgzg′z∂2uj+1∂x̂2−∂2uj+1∂ẑ2−γy1+γgzg′z∂2uj+1∂ỹ2−∂2uj+1∂z̃2+γ2xy1+γgz2∂2uj+1∂x̄2−∂2uj+1∂ȳ2−g″zg′z3∂uj+1∂z+1g′z2∂2uj+1∂z2−γx1+γgz∂ψzj+1∂x−γy1+γgz∂ψzj+1∂y+1g′z∂ψzj+1∂z],(10g)$

where uj = u(x, y, z, tj) represents the scalar wavefield at the jth time step in the trapezoid coordinate system; (xs, ys, zs) is the position of the source in the trapezoid coordinate system;

$στ=3vmax2Lττ̄Lτ2$

ln

$1R$

, where Lτ is the thickness of CPML absorbing boundary along the τ direction,

$τ̄$

denotes the distance to the inner area in the τ direction; and

$ατ=αmax1−τ̄Lτ$

, τ ∈ {x, y, z}.

A schematic of the grid discretization in the 3D trapezoid-grid CPML area is shown in Figure 3. In this work, 30 and 20 absorbing boundary layers are usually used for the trapezoid-grid CPML area in the horizontal and vertical directions, respectively.

FIGURE 3. Schematic of the grid discretization in the 3D trapezoid-grid CPML area.

### 2.3 Stability Analysis

Stability condition is usually required for the FD scheme to give a stable time step. From Eq. 10a, we use a local frozen coefficients technique in each discrete point and can get the full-discretization form of the 3D trapezoid coordinate system acoustic equation without the CPML boundary condition and source function:

$1v2um,n,lj+1−2um,n,lj+um,n,lj−1Δt2−γ2xm2+11+γgzl21Δ2∑p=1Nxηpxum−p,n,lj−2um,n,lj+um+p,n,lj−γ2yn2+11+γgzl21Δ2∑p=1Nyηpyum,n−p,lj−2um,n,lj+um,n+p,lj−2γ2xm1+γgzl21Δ∑p=1Nxcpxum+p,n,lj−um−p,n,lj−2γ2yn1+γgzl21Δ∑p=1Nycpyum,n+p,lj−um,n−p,lj+γxm1+γgzlg′zl12Δ2∑p=1Nxzηpxzum−p,n,l−pj+um+p,n,l+pj−um−p,n,l+pj−um+p,n,l−pj+γyn1+γgzlg′zl12Δ2∑p=1Nyzηpyzum,n−p,l−pj+um,n+p,l+pj−um,n−p,l+pj−um,n+p,l−pj−γ2xmyn1+γgzl212Δ2∑p=1Nxyηpxyum−p,n−p,lj+um+p,n+p,lj−um−p,n+p,lj−um+p,n−p,lj+g″zlg′zl31Δ∑p=1Nzcpzum,n,l+pj−um,n,l−pj−1g′zl21Δ2∑p=1Nzηpzum,n,l−pj−2um,n,lj+um,n,l+pj=0,(11)$

where

$um,n,lj$

is the wavefield at (xm, yn, zl, tj), xm = x0 + (m − 1)Δx, yn = y0 + (n − 1)Δy, zl = z0 + (l − 1)Δz, tj = t0 + (j − 1)Δt, Nx, Ny, Nz, Nxy, Nxz, Nyz are half-of-spatial FD orders, ηx, ηy, ηz, ηxy, ηxz, ηyz are corresponding FD coefficients of the second-order derivative, and cx, cy, cz, cxy, cxz, cyz are corresponding FD coefficients of the first-order derivative.

To derive the stability condition, we use the plane wave solution that is defined as

$ux,y,z,t=u0*eiωt−ikxx−ikyy−ikzz,(12)$

where

$u0*$

is the amplitude of the plane wave, i is the imaginary unit, ω is the angular frequency, and kx, ky, kz are wavenumbers in the x-, y- and z-directions, respectively. Similar to stability analysis of Kosloff and Baysal (1982), by substituting Eq. 12 into Eq. 11 and only considering the maximum wavenumber, the stability condition of the 3D acoustic equation in the trapezoid coordinate system can be expressed as

$Δt<Δvmaxmaxm,n,lγ2xm2+11+γgzl2∑p=1Nxmodp,2ηpx+γ2yn2+11+γgzl2∑p=1Nymodp,2ηpy+1g′zl2∑p=1Nzmodp,2ηpz−12,(13)$

where mod is the function for the getting remainder, and

$maxm,n,l$

represents the maximum value of the objective function at those discrete points (xm, yn, zl).

## 3 Numerical Results

In the following numerical examples, Eq. 10 is discretized by the eighth-order FD in the space, and conventional Taylor-expansion–based high-order FD coefficients (Dablain, 1986) are adopted.

### 3.1 Homogenous Model

First, we use a 3D homogenous model with a constant velocity of 2000 m/s to verify the effectiveness of our trapezoid-grid FDTD method and corresponding CPML absorbing boundary condition. A Ricker wavelet with a dominant frequency of 20 Hz is located at the center of the model as the source. The FD time step is taken as 1.6 ms. The scaling parameter γ is set as 2.78 × 10–4, the sampling function g(z) = z, and the lateral grid sizes in the Cartesian coordinate system increase from 7.5 to 10 from top to bottom. Figure 4A shows the snapshot obtained by the trapezoid-grid FDTD with CPML at 0.45 s, while Figure 4B shows the corresponding snapshot without CPML. Figure 5 shows the comparison between the recorded seismograms computed by the uniform-grid FDTD and the trapezoid-grid FDTD at (x0, y0, z0) = (600, 600, 0 m). The comparison between Figures 4A,B demonstrates that trapezoid-grid CPML can effectively reduce boundary reflections, while Figure 5 demonstrates the accuracy of the trapezoid-grid FDTD method for the homogenous model.

FIGURE 4. Snapshots obtained by the trapezoid-grid FDTD for the homogenous model (A) with CPML and (B) without CPML.

FIGURE 5. Comparison of seismograms obtained by the uniform-grid FDTD and the trapezoid-grid FDTD for the homogenous model at (x0, y0, z0) = (600 m, 600 m, 0 m). The black solid line represents the uniform-grid FDTD result, and the red dash line represents the trapezoid-grid FDTD result.

### 3.2 Overthrust Model

Then, we apply our method to the SEG/EAGE overthrust model (Figure 6A), which is based on the real overthrusts of South America. Figures 6B,C show the modeling area of the SEG/EAGE overthrust model in the trapezoid coordinate system and the Cartesian coordinate system, respectively. A Ricker wavelet with the dominate frequency of 4.2 Hz is located at (10 km, 10 km, 0.5 km) as the source. The grid sizes for the uniform-grid FDTD are 50 m × 50 m × 50 m, which means the minimum NPPW in each direction is close to 10. We therefore set the minimum NPPW in x0-, y0-, and z0-direction as 10 for the trapezoid-grid FDTD method, and get the scaling parameter as γ = 2.07 × 10–4. Figure 7A shows the vertical sampling function g(z) used for this model, and the vertical grid sizes in the Cartesian coordinate system increase from 51.8 m in the shallow region to a maximum value of 142.6 m in the deep region, as shown in Figure 7B. Based on the stability analysis, the time step for the trapezoid-grid FDTD and the uniform-grid FDTD are 3.697 and 3.585 ms, respectively. Receivers are located on the surface along y0 = 10 km. Figure 8A shows the snapshot at 2.5 s computed by the uniform-grid FDTD, and Figure 8B is the snapshot at 2.5 s in the trapezoid coordinate system computed by our trapezoid-grid method. Using coordinate transformation and cubic spline interpolation, we can get the corresponding snapshot in the Cartesian coordinate system, as shown in Figure 8C. Figures 8A,C show good agreement. To give more detailed comparisons, single-trace seismograms at (7.5, 10, 0 km), (10, 10, 0 km), and (12.5, 10, 0 km) for both the uniform-grid (black solid line) and the trapezoid-grid (red dash line) FDTD are shown in Figure 9. Figure 9 also shows good agreement between the uniform-gird FDTD and the trapezoid-grid FDTD. On our computing platform (Intel(R) Xeon(R) Sliver 4216 CPU @ 2.10GHz, 256GB of memory, and C++ codes), using 16-threads computation and similar code optimization techniques, the running time for the trapezoid-grid FDTD and the uniform-grid FDTD is calculated as 2203 s and 2925 s, respectively, which shows a calculation efficiency improvement of 24.7%. The memories for the trapezoid-grid FDTD and the uniform-grid FDTD are about 336 and 1213 MB, respectively, which shows a memory reduction of 72.3%. Considering that the simulation area of our trapezoid-grid method is almost 60% of that of the uniform-grid method, for the common simulation area, we can achieve about 50% reduction on memory usage.

FIGURE 6. (A) SEG/EAGE overthrust model; (B) SEG/EAGE overthrust model in the trapezoid coordinate system; (C) actual simulating area for the trapezoid-grid FDTD in the Cartesian coordinate system.

FIGURE 7. (A) Vertical sampling function g(z); (B) variation of Δz0 for the trapezoid-grid FDTD.

FIGURE 8. Snapshots for the SEG/EAGE overthrust model at 2.5 s: (A) uniform-grid FDTD result; (B) trapezoid-grid FDTD result in the trapezoid coordinate system; (C) trapezoid-grid FDTD result in the Cartesian coordinate system.

FIGURE 9. Single-trace seismograms for the SEG/EAGE overthrust model at the receiver position of (A) (7.5, 10, 0 km); (B) (10, 10, 0 km); and (C) (12.5, 10, 0 km). The black solid lines represent the uniform-grid FDTD results, and the red dash lines represent the trapezoid-grid FDTD results.

## Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

## Author Contributions

BW contributed to the conception of the study. Writing, compiling, and debugging of the program were completed by WX. WT wrote the first draft of the manuscript, and all authors contributed to the revision.

## Funding

We would like to thank the Natural Science Foundation of Shaanxi Province under Grant No. 2020JM-018 and National Natural Science Foundation (41974122) of China for funding this work.

## Conflict of Interest

Authors BW, WT, and BL are employed by SINOPEC.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.