## 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 2*M* 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

In this article, the 3D trapezoid coordinate transformation is defined as

where (*x*_{0}, *y*_{0}, *z*_{0}) 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 z_{0}-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:

where *f*_{0} is the dominant frequency of the source term, *N*_{0} is the preferred NPPW and is related to the accuracy of the adopted FD scheme, and *v*_{min}(*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 *v*_{min}(*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 *v*_{min}(*g*(*z*)). In particular, we usually set the order of the local polynomial as three. Such *v*_{min}(*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 z_{0}-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 Δ*x*_{0}(*z*) and Δ*y*_{0}(*z*) are always smaller than or equal to Δ*z*_{0}(*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

where *u*^{j} = *u*(*x*_{0}, *y*_{0}, *z*_{0}, *t*_{j}) represents the scalar wavefield at the *j*th time step in the Cartesian coordinate system; *v* is the velocity; Δ*t* is the time interval; *f*(*t*) is the source term;

is the position of source;

ln

, where R denotes the designated theoretical boundary reflection coefficient, *v*^{max} is the maximum velocity of the model,

is the thickness of CPML absorbing boundary along the *τ*_{0} direction, and

denotes the distance to the inner area in the *τ*_{0} direction;

and *α*^{max} = *πf*_{0}.

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

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

where

and

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

and

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

and

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

axes, *θ*_{2} is the angle between *y* and

axes, and *θ*_{3} is the angle between *x* and

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

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

. 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

where *u*^{j} = *u*(*x*, *y*, *z*, *t*_{j}) represents the scalar wavefield at the *j*th time step in the trapezoid coordinate system; (*x*^{s}, *y*^{s}, *z*^{s}) is the position of the source in the trapezoid coordinate system;

ln

, where *L*_{τ} is the thickness of CPML absorbing boundary along the *τ* direction,

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

, *τ* ∈ {*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.

### 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:

where

is the wavefield at (*x*_{m}, *y*_{n}, *z*_{l}, *t*_{j}), *x*_{m} = *x*^{0} + (*m* − 1)Δ*x*, *y*_{n} = *y*^{0} + (*n* − 1)Δ*y*, *z*_{l} = *z*^{0} + (*l* − 1)Δ*z*, *t*_{j} = *t*^{0} + (*j* − 1)Δ*t*, *N*_{x}, *N*_{y}, *N*_{z}, *N*_{xy}, *N*_{xz}, *N*_{yz} are half-of-spatial FD orders, *η*^{x}, *η*^{y}, *η*^{z}, *η*^{xy}, *η*^{xz}, *η*^{yz} are corresponding FD coefficients of the second-order derivative, and *c*^{x}, *c*^{y}, *c*^{z}, *c*^{xy}, *c*^{xz}, *c*^{yz} are corresponding FD coefficients of the first-order derivative.

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

where

is the amplitude of the plane wave, *i* is the imaginary unit, *ω* is the angular frequency, and *k*_{x}, *k*_{y}, *k*_{z} 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

where *mod* is the function for the getting remainder, and

represents the maximum value of the objective function at those discrete points (*x*_{m}, *y*_{n}, *z*_{l}).

## 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 (x_{0}, y_{0}, z_{0}) = (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 (x_{0}, y_{0}, z_{0}) = (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 x_{0}-, y_{0}-, and z_{0}-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 y_{0} = 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 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.

## 4 Conclusion

In this article, we propose a 3D trapezoid-grid FDTD seismic wave modeling method based on the increasing trend of seismic wave velocity with depth. The trapezoid-grid mesh in the physical Cartesian system can effectively reduce the oversampling in the high-velocity region compared with the uniform-grid method, and the design of 3D trapezoid coordinate transform greatly avoids the difficulty of processing an irregular grid. We derive the 3D acoustic equation in the trapezoid coordinate system. The corresponding CPML boundary condition is also given to decrease artificial boundary reflection. To obtain a stable and efficient wave modeling result, we combine the plane wave theory and frozen coefficients technique and provide an effective stability condition for the 3D trapezoid-grid FDTD method. The discretization of the 3D acoustic equation in the trapezoid coordinate system is completed by the eighth order and second order finite-difference method in the space and time domain, respectively. The 3D homogenous model is given to verify the effectiveness of trapezoid-grid FDTD and the performance of the CPML boundary. Numerical tests on the SEG/EAGE overthrust model indicate the accuracy and the significant memory reduction of our method compared with uniform-grid FDTD. The key idea of our method is the combination of the trapezoid coordinate transformation and the FD stencils. Such idea can be generalized to many other wave equations such as elastic equation (Zhan et al., 2017) and Maxwell’s equations (Zhan et al., 2021). Besides, our method is actually dealing with the regular grids in the trapezoid coordinate system, which means that we can combine other methods to treat the irregular surface (Li et al., 2020) or curved interfaces (Zhan et al., 2020).

## 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.

## References

Abreu, R., Stich, D., and Morales, J. (2015). The Complex-Step-Finite-Difference Method. *Geophys. J. Int.* 202, 72–93. doi:10.1093/gji/ggv125

##### Disclaimer:

This article is autogenerated using RSS feeds and has not been created or edited by OA JF.

Click here for Source link (https://www.frontiersin.org/)