https://orcid.org/0000-0002-9341-520X
SUMMARY
The finite-difference method (FDM), limited by uniform grids, often encounters severe oversampling in high-velocity regions when applied to multiscale subsurface structures, leading to reduced computational efficiency. A feasible solution to this issue is the use of non-uniform grids. However, previous discontinuous grid approaches required careful consideration of interpolation operations in transition regions, while single-block continuous grids lacked flexibility. This paper proposes a novel approach using multiblock stretched grids with positive and negative singularities to achieve non-uniform grids, the numerical simulation of seismic waves is realized by combining it with the curvilinear grid FDM (CGFDM). Our method facilitates seamless information exchange between coarse and fine grids without additional interpolation or data processing and allows for flexible grid configurations by adjusting singularity pairs.
The effectiveness of our approach is verified through comparisons with the generalized reflection/transmission method and the finite-element method. Numerical experiments demonstrate the method's reliable accuracy and significant reduction in grid points compared to uniform grids. Although the stability of our method has not been rigorously mathematically proven, we demonstrate that the algorithm remains applicable for sufficiently long simulations to address realistic scenarios.
INTRODUCTION
In the real Earth, complex medium characteristics, such as heterogeneity, anisotropy and irregular boundaries, present significant challenges for seismic wave numerical simulations. These features lead to diverse seismic wave propagation behaviours, increasing the complexity of simulations. While the finite-difference method (FDM) is widely favoured for its simplicity and computational efficiency (Graves 1996; Moczo et al. 2007), it is generally applied to uniform grids, which proves insufficiency when addressing multiscale underground structures, such as low-velocity layers (Zorin et al. 2002) and fractured-cavernous reservoirs (Yang et al. 2012). Sufficiently accurate simulations require grid refinement based on minimum velocity, which can increase memory usage and reduce efficiency. A practical solution is to use a fine grid in targeted zones while employing a coarse grid in less critical areas (Moczo 1989; Chu & Stoffa 2012). Thus, exploring flexible grid partitioning methods is essential. Non-uniform grids allow for refinement based on geological complexity, optimizing grid distribution and enhancing computational efficiency.
Masson (2023) proposed a distributional finite-difference method (DFDM) based on the weak form of partial differential equations (PDEs), enabling the use of non-conformal grids for non-uniform grid configurations (Masson & Virieux 2023; Masson et al., 2024). Here, we primarily focus on the more straightforward approach of the FDM based on the strong form of PDEs. The FDM offers two primary means for implementing non-uniform grids, one of which involves stable discontinuous grids, as illustrated in Fig. 1(a). This technique, introduced by Kristek et al. (2010), utilizes various interpolation methods to facilitate information exchange between coarse and fine grids (Zhang et al. 2013; Li et al. 2015), allowing for grid refinement and significantly enhancing computational efficiency. Simultaneously, to ensure long-term stability and accuracy, this approach incorporates filtering operations. Dynamic injection (Miao & Zhang 2023) and varying weighting coefficients (Huang & Dong 2009) using Taylor expansion can facilitate information exchange between grids. However, Taylor expansion achieves high accuracy only when variables are near the expansion point (Zhang & Yao 2017), which may lead to numerical dispersions in the FDM. In addition, another method involves cleverly setting up the grid and using rotated finite-difference operators in transition regions to connect coarse and fine grids (Moczo et al. 1996). Yet, this method might present challenges in maintaining long-term stability (Fan et al. 2015). The second method involves the continuous grid method, as shown in Fig. 1(b), which entails uniformly stretching the grid in a specific direction to achieve non-uniform grid configurations (Zhang et al. 2012a; Li & Jia 2018). Continuous grid method does not require interpolation or similar operations for information transfer, providing higher accuracy and stability. However, this overall stretching can lead to severe distortions in certain areas of the grid. Therefore, we need other, more effective strategies to achieve flexible non-uniform grids.
Two types of non-uniform grids: (a) discontinuous grid, (b) continuous grid.
Previous continuous grid methods were typically based on stretching individual grid blocks. Moczo (1989) combined a FDM using a rectangular grid with non-uniform grid spacing to compute 2-D SH waves. This method achieved higher efficiency compared to conventional approaches in media where internal velocity and density vary linearly in the horizontal and vertical directions. Building on this concept, Zhang et al. (2012a) proposed a curvilinear grid finite-difference method (CGFDM), which used non-uniform grids stretched along the depth direction to solve seismic wave equations in spherical coordinates, improving the computational efficiency of the FDM. However, this approach could not stretch the grid in other directions. Pitarka (1999) introduced a method akin to Moczo (1989), employing a fourth-order staggered-grid finite-difference operator for simulating 3-D elastic waves on non-uniform grids, thereby enhancing the efficiency of finite fault simulations. Sofronov et al. (2015) also conducted similar work, using non-uniform grids in two directions for elastic dynamic simulations with a multiblock FDM. Nevertheless, this approach resulted in excessive stretching in non-target regions. In addition, trapezoidal grid FDM (Wu et al. 2019, 2022) is another type of non-uniform grid. Trapezoidal grids allow the grid spacing to increase with depth, aligning with the tendency for seismic wave velocity to increase with depth. This reduces oversampling in high-velocity regions. Nevertheless, as the grid spacing increases with depth, the horizontal spacing in the simulation domain also increases. These grid-stretching techniques based on individual grid blocks have limitations when dealing with complex Earth media. Thus, we need to explore more flexible methods, such as multiblock grids stretching, to overcome these challenges more effectively.
In the early stages of computational fluid dynamics (CFD), the concept of multiblock grids was introduced to overcome the limitations of single-block grids in addressing fluid problems with complex boundaries (Enwald et al. 1999; Badcock et al. 2000; Ruiz-Gironés & Sarrate 2010). By dividing the computational domain into multiple subregions, this approach allows for more flexible and accurate simulations of intricate geometries and flow patterns (Yu et al. 2002). Additionally, multiblock structured grids offer several advantages in numerical analysis, including decomposing the grid into blocks with simpler topological structures, enabling direct numerical analysis of complex geometries and flows (Blazek 2015). Each block can be solved independently, making the multiblock method particularly suitable for parallel computing, significantly reducing computational time (Takaki et al. 2003).
We present a type of multiblock structured grid (Weatherill & Forsey 1985; Wang et al. 1998) where the nodes between adjacent blocks are aligned, ensuring seamless information exchange between neighbouring grid blocks without complex interpolation or data processing. This method is straightforward to implement, promoting efficient communication, maintaining computational consistency across grid blocks, and enhancing overall simulation efficiency. The computational domain is divided into four-sided blocks, with adjacent blocks connected at nodes of regular connectivity, except at singular points. An interior node (Sun et al. 2021) is classified as a negative singularity if it is bounded by three quadrilateral grid cells, and as a positive singularity if it is bounded by five quadrilateral grid cells, as shown in Fig. 2. These two types of singularities are most commonly used in grid generation techniques, and therefore, this paper focuses exclusively on these types. Fogg et al. (2018) introduced a method for controlling grid density variation in multiblock grids by adding pairs of positive and negative singularities, a technique we propose applying in the finite-difference simulation of seismic waves. By introducing these singularities, we can adaptively refine the grid (Fogg et al. 2015; Armstrong et al. 2019) conformly with the complexity of seismic wave propagation. This approach aims to enhance both computational efficiency and accuracy, addressing the challenges posed by diverse subsurface structures in seismic wave simulations. When dealing with singularities in the grid, a finite-difference operator based on summation-by-parts (SBP) (Fernández et al. 2014; Svärd & Nordström 2014) is often combined with simultaneous approximation terms (SAT) techniques (Sun et al. 2020) to facilitate information exchange between multiple grid blocks. The SBP-SAT method lies in its ability to establish strict stability proofs using the energy method (Dovgilovich & Sofronov 2015), leading to robustness in spatial and time domains. However, the application of SAT techniques can reduce boundary accuracy and increase implementation complexity. In this work, we use the curvilinear grid finite-difference method (CGFDM) proposed by Zhang & Chen (2006) and Zhang et al. (2012b), which simplifies the process by allowing direct information exchange through the matching of grid points, avoiding the need for complex interpolation techniques.
Controlling grid density with an extra singularity pair and examples of grid singularities in a multiblock decomposition.
This paper presents a novel non-uniform grid FDM integrating multiblock grids with singularity pairs and the CGFDM. First, we briefly introduce the grid generation techniques required for this approach. Then, we review the numerical implementation of the CGFDM for solving elastic wave equations. Following that, we explain how CGFDM handles positive and negative singularities. Next, we apply our proposed method to simulate seismic wave propagation on non-uniform grids. Finally, we conduct a series of numerical experiments to verify the effectiveness of our non-uniform grid FDM.
METHODOLOGY
Grid generation
In 2-D, several methods have been developed for generating quadrilateral (quad) grids, including paving (Blacker & Stephenson 1991), Cartesian grid (Schneiders 1996), submapping (White 1996) and medial-axis-based decomposition methods (Rigby 2003). These algorithms can automatically generate multiblock topological structures for complex geometries. However, the geometries required for simulating seismic wave propagation in this research are not overly complex. As a result, manually constructing multiblock structured grids is a common approach, as it creates higher quality grids and provides more precise control over their shape and topology.
In industry, the standard method for creating multiblock grids for complex shapes involves manually generating the blocks using tools such as ICEM or Pointwise. This approach is typically carried out in two steps. First, the positions of singular points are identified, after which blocks are generated based on these singular points, with the block edges and vertices aligned to the actual geometry. As shown in Fig. 2, a pair of positive and negative singularities are arranged, and the grid blocks are divided based on these singular points, resulting in nine subgrid blocks. The grid density gradually increases from the positive singularities towards the negative singularities.
At the boundaries of the grid in Fig. 2, good orthogonality is maintained, which facilitates the application of perfectly matched layer (PML) boundary conditions. Furthermore, multiple grids of this kind can be assembled to create more flexible non-uniform grid configurations. In this study, singularity pairs are used to enable non-uniform grid generation. For more complex geometries, alternative singularity configurations (Armstrong et al. 2015; Fogg et al. 2015) could be applied, though such cases are beyond the scope of this work.
Velocity–stress equation and curvilinear grid finite-difference method
In 2-D Cartesian coordinates, the velocity–stress equations for elastic wave propagation can be written in a compact form using vector and tensor notation:
where |${{\bf v}} = {{[ {{{v}_x},{{v}_y}} ]}^T}$| is the velocity vector,
We apply the CGFDM to solve the aforementioned elastic wave equations in the curvilinear coordinate system for each grid block, which can be considered a boundary-defined body-fitted grid. This method is particularly effective in handling complex geometries, as it transforms the physical domain into a computational domain that better aligns with the grid structure, ensuring a more accurate representation of the wave equations even in domains with irregular boundaries. This scheme has been successfully applied to simulations of strong ground motion (Zhang et al. 2012b) and earthquake rupture dynamics (Zhang et al. 2014). For more detailed information, refer to Zhang & Chen (2006) and Zhang et al. (2012b).
The CGFDM uses a collocated grid with the MacCormack finite-difference operator and fourth-order Runge–Kutta time integration. The MacCormack-type scheme alternates between forward and backward one-sided differences during the multistage time integration. In this study, we employ the DRP/opt (Dispersion Relation Preserving/optimized) MacCormack operator (Hixon 1997), which optimizes the scheme's dispersion and dissipation errors, significantly improving grid resolution. The spatial derivatives with respect to ξ-axis can be approximated by the forward and backward finite-difference operators:
where Lξ represents the spatial difference with respect to the ξ-axis; the subscript i is the grid index; the superscript F and B denote the forward and backward operators, respectively; and a−1 = −0.30874, a0 = −0.6326, a1 = 1.2330, a2 = −0.3334 and a3 = 0.04168 (Hixon 1997).
Additionally, two types of boundary conditions are implemented: a free-surface boundary condition at the top, achieved through a traction image method (Zhang & Chen 2006), and an absorbing boundary condition using the PML (Berenger 1994), specifically the unsplit complex frequency-shifted PML (Zhang & Shen 2010), which enhances wave absorption at the boundaries.
Information exchange between grid blocks and computation at singularity points
In the case of Fig. 2, CGFDM can be easily applied within each grid block, and since the grid points between adjacent blocks are aligned (except for a few points near singularities), point-to-point information exchange between neighbouring grid blocks is facilitated using ghost points.
As shown in Fig. 3, handling the computation around positive and negative singularities and their surrounding grid points is remarkably straightforward. For a positive singularity, located at the intersection of five grid lines, one of the grid lines is removed, leaving four lines to construct the difference stencil in both computational directions. For a negative singularity, located at the intersection of three grid lines, one of the grid lines is reused, allowing four grid lines to be employed for stencil construction. Since the DRP/opt MacCormack operator has a stencil length of 7, some points near the singularity will require the singularity's value to complete their difference stencils. The grid points needed for stencil construction at these points are explained in detail in the tables of Figs 3(d) and (e). In these tables, ‘r’ indicates that the stencil is completed along the red grid lines, while ‘b’ indicates completion along the blue grid lines.
Illustration of the computation at singularity points. (a) Grid lines near a positive singularity. (b) The difference stencil in the computational domain, with a stencil length of 7. (c) Grid lines near a negative singularity. In (a) and (c), the blue grid lines indicate that the difference stencil requires the use of the singularity, while the red grid lines represent the directions where the stencil is applied normally. (d) and (e) show the specific points needed to construct the difference stencils at positive and negative singularities, respectively. The tables also illustrate part of the difference stencils that require the use of the singularity grid points, where ‘r’ represents completing the stencil along the red grid lines, and ‘b’ represents completing the stencil along the blue grid lines.
NUMERICAL EXPERIMENT
In this section, three numerical experiments are conducted to demonstrate the effectiveness of our method. The reference solutions used in these experiments include the generalized reflection/transmission coefficient method (GRTM; Luco & Apsel 1983; Chen 1993) and the finite-element method (FEM). The GRTM, by solving wave equations in a semi-analytical framework, provides analytical reference solutions for homogeneous full-space and layered models, while FEM is used for more complex geometries. The FEM reference solutions are generated using the commercial software COMSOL Multiphysics, which employs a time-domain explicit FEM based on third-order Lagrange polynomials. Spatial discretization is achieved with unstructured quadrilateral meshes, and time integration is performed using a fourth-order Runge–Kutta scheme. COMSOL Multiphysics is renowned for its user-friendly interface, making it accessible and efficient for complex simulations.
The accuracy of the solutions is evaluated through global measurement, including envelope misfit (EM), phase misfit (PM), envelope goodness-of-fit (EG) and phase goodness-of-fit (PG), based on the time-frequency (TF) misfit criteria proposed by Kristeková et al. (2006; 2009). EM and PM quantify the discrepancies between the solutions in amplitude and phase, respectively, while EG and PG measure the agreements in these terms. These TF criteria effectively differentiate between amplitude and phase errors and are widely used in seismological assessments of numerical solutions (Fichtner 2010; Nissen-Meyer et al. 2014; Sun et al. 2017). According to Kristeková et al. (2009), the ‘excellent’ thresholds for EM and PM are 0.22 and 0.2, respectively, meaning EM or PM values should fall below these limits. For EG and PG, the minimum ‘excellent’ level is 8, so these values should exceed this threshold.
For the experiments, a Ricker wavelet is used as the source time function, expressed as:
where fc is the central frequency, and t0 is the time delay. The points per wavelength (PPW) are calculated as vs/(2.5fchs), where vs is the minimum velocity, and hs is the maximum spatial interval between adjacent grid nodes in both x- and y-directions. The source is implemented as a vertical point force applied to a single grid node. The results are directly compared with the reference solutions without applying any filtering or post-processing operations.
We use the effective media parameters approach (Moczo et al. 2002) to account for internal discontinuous interfaces. For a more accurate consideration of interface issues, one can refer to recent updates on the FDM (Moczo et al. 2023; Valovcan et al. 2023, 2024), which provide solutions for more precisely handling interfaces. These updates delve into the impact of the wavenumber bandwidth limitation of the medium and the corresponding spatial sampling constraints of the grid on the accurate representation of interfaces in numerical simulations.
Homogeneous full space model
In the first numerical experiment, we consider a homogeneous full-space model. Two grid blocks are constructed as shown in Fig. 2, which are combined to form the grid depicted in Fig. 4. This grid contains 537 182 grid points, representing a 35.14 per cent reduction compared to the 828 259 grid points required for a uniformly fine grid. We perform three simulations differing from each other by position of the source with parameters fc = 2.50 Hz and t0 = 0.5657 s: the fine grid (S1, x : 4040.302 m, y : −1323.60 m), the transition grid area (S2, x : 4238.78 m, y: −3979.27 m) and the coarse grid (S3, x : 4116.84 m, y : −6291.51 m). The medium is homogeneous elastic with a P-wave velocity of 2.0 km s−1, an S-wave velocity of 1.0 km s−1 and a density of 1.0 g cm−3. The coarse grid spacing is 9.12 m and the fine grid spacing is 4.56 m, yielding a PPW of 11.91. The time-step for the simulation is set at 0.0015 s, with a total of 5 000 time-steps. The grid has orthogonal boundaries, and 12 layers of PML boundary conditions are applied around the perimeter. Four receivers are sparsely distributed among the fine, transition and coarse grids, with their specific locations referenced in the ‘Experiment 1’ column of Table 1. Notably, Receivers 2 and 3 are located near the singularity.
Snapshots of the wavefields for vx and vy at t = 3.00 s when the source is located in the finer grid region. The blue grid lines represent the computational grid used for the simulation, and the grid density has been downsampled for clarity. The black five-pointed star indicates the source location, and the multiple black triangles mark the positions of different receivers. Subsequent figures follow similar notation and representation.
Detailed coordinates of all receivers in the three numerical experiments presented in this paper (Units: m).
| No. | Experiment 1 (x, y) | Experiment 2 (x, y) | Experiment 3 (x, y) |
|---|---|---|---|
| Receiver 1 | (2333.55, −1433.55) | (4950.43, 0.00) | (3997.27, 0.00) |
| Receiver 2 | (2326.79, −2863.26) | (8133.62, 0.00) | (4997.27, 0.00) |
| Receiver 3 | (2397.83, −4339.99) | (11015.09, 0.00) | (5997.79, 0.00) |
| Receiver 4 | (2356.45, −6079.58) | (4946.58, −4433.97) | (6999.36, 0.00) |
| Receiver 5 | \ | (8144.10, −4469.10) | (3996.15, −503.77) |
| Receiver 6 | \ | (11019.69, −4459.54) | (5008.43, −502.47) |
| Receiver 7 | \ | (4951.32, −7545.67) | (5993.47, −510.25) |
| Receiver 8 | \ | (8152.98, −7544.01) | (6999.25, −504.11) |
| Receiver 9 | \ | (11000.17, −7544.51) | \ |
| No. | Experiment 1 (x, y) | Experiment 2 (x, y) | Experiment 3 (x, y) |
|---|---|---|---|
| Receiver 1 | (2333.55, −1433.55) | (4950.43, 0.00) | (3997.27, 0.00) |
| Receiver 2 | (2326.79, −2863.26) | (8133.62, 0.00) | (4997.27, 0.00) |
| Receiver 3 | (2397.83, −4339.99) | (11015.09, 0.00) | (5997.79, 0.00) |
| Receiver 4 | (2356.45, −6079.58) | (4946.58, −4433.97) | (6999.36, 0.00) |
| Receiver 5 | \ | (8144.10, −4469.10) | (3996.15, −503.77) |
| Receiver 6 | \ | (11019.69, −4459.54) | (5008.43, −502.47) |
| Receiver 7 | \ | (4951.32, −7545.67) | (5993.47, −510.25) |
| Receiver 8 | \ | (8152.98, −7544.01) | (6999.25, −504.11) |
| Receiver 9 | \ | (11000.17, −7544.51) | \ |
Detailed coordinates of all receivers in the three numerical experiments presented in this paper (Units: m).
| No. | Experiment 1 (x, y) | Experiment 2 (x, y) | Experiment 3 (x, y) |
|---|---|---|---|
| Receiver 1 | (2333.55, −1433.55) | (4950.43, 0.00) | (3997.27, 0.00) |
| Receiver 2 | (2326.79, −2863.26) | (8133.62, 0.00) | (4997.27, 0.00) |
| Receiver 3 | (2397.83, −4339.99) | (11015.09, 0.00) | (5997.79, 0.00) |
| Receiver 4 | (2356.45, −6079.58) | (4946.58, −4433.97) | (6999.36, 0.00) |
| Receiver 5 | \ | (8144.10, −4469.10) | (3996.15, −503.77) |
| Receiver 6 | \ | (11019.69, −4459.54) | (5008.43, −502.47) |
| Receiver 7 | \ | (4951.32, −7545.67) | (5993.47, −510.25) |
| Receiver 8 | \ | (8152.98, −7544.01) | (6999.25, −504.11) |
| Receiver 9 | \ | (11000.17, −7544.51) | \ |
| No. | Experiment 1 (x, y) | Experiment 2 (x, y) | Experiment 3 (x, y) |
|---|---|---|---|
| Receiver 1 | (2333.55, −1433.55) | (4950.43, 0.00) | (3997.27, 0.00) |
| Receiver 2 | (2326.79, −2863.26) | (8133.62, 0.00) | (4997.27, 0.00) |
| Receiver 3 | (2397.83, −4339.99) | (11015.09, 0.00) | (5997.79, 0.00) |
| Receiver 4 | (2356.45, −6079.58) | (4946.58, −4433.97) | (6999.36, 0.00) |
| Receiver 5 | \ | (8144.10, −4469.10) | (3996.15, −503.77) |
| Receiver 6 | \ | (11019.69, −4459.54) | (5008.43, −502.47) |
| Receiver 7 | \ | (4951.32, −7545.67) | (5993.47, −510.25) |
| Receiver 8 | \ | (8152.98, −7544.01) | (6999.25, −504.11) |
| Receiver 9 | \ | (11000.17, −7544.51) | \ |
Wavefield snapshots at t = 3.00 s are shown in Fig. 4, Fig. 5 and Fig. 6, where P and S waves excited by the point source are clearly visible. These waves propagate stably across the fine, transition and coarse grids. Fig. 7 compares waveforms recorded at four receivers for different source locations against the GRTM solution for a homogeneous full space. The waveforms from both methods align closely, with EM and PM values remaining below 1 per cent, and EG and PG values exceeding 9.9. Based on time-frequency misfit analysis, all metrics fall within the ‘excellent’ category, indicating that our method produces results that are in strong agreement with the reference solution. We extended the number of time steps to 600 000 for the case where the source is located in the fine grid (S1), and Fig. 8 shows the waveform recorded at Receiver 2. Our method remains stable even after such a long simulation time, with no signs of exponential growth. To further test the stability of our approach, we swapped the positions of Source S1 and Receiver 4, adjusting the direction of the source to be opposite to its original direction to ensure waveform consistency. Fig. 9 shows the comparison results of the waveforms simulated by both, with high matching accuracy and an error margin of approximately one per cent, indicating that our method satisfies simple numerical reciprocity.
Snapshots of the wavefields for vx and vy at t = 3.00 s when the source is located in the transition grid region.
Snapshots of the wavefields for vx and vy at t = 3.00 s when the source is located in the coarser grid region.
Synthetic seismograms of the vx and vy components at the receivers shown in Figs 4, 5 and 6. The red solid lines represent the GRTM solution, while the blue dashed lines correspond to the CGFDM solution. The left-side labels indicate the receiver numbers, and the misfit quantification for the CGFDM is marked near each waveform. The same labelling convention applies to subsequent figures.
Waveforms of vx and vy at Receiver 2 from Fig. 4, in a homogeneous full-space model, with the source located at S1. The simulation was run for a total of 600 000 time-steps, corresponding to a total simulation time of 900 s.
Comparison of waveforms simulated before and after swapping the positions of the source S1 and Receiver 4. The solid red line represents the waveform simulated when the source is at S1 and the receiver at Receiver 4, while the dashed blue line shows the opposite scenario.
To assess the stability of the multiblock grid configuration, we also examined wavefield propagation with the source located at singularity points. Fig. 10 illustrates the energy evolution over time across the entire simulation region when the source is located at singularities, specifically at Receiver 2 and Receiver 3 as shown in Fig. 4. The main energy of the wavefield is attenuated after approximately 6 000 time-steps, and our method ensures that the energy in the entire simulation region does not increase significantly before 600 000 steps, which meets the simulation requirements of most practical problems. The results show that even when the source is placed at a singularity, our method remains stable within a certain scale. In addition to these cases, we tested several other source locations, and the results remained stable. However, when the source is positioned too close to the singularity, the deviation from the reference values increases. Although configuring the multiblock grids may introduce some instability, these numerical tests demonstrate that our method can maintain stability over realistic simulation durations.
The evolutions of energy over time in the entire simulation domain when the source is located at positive and negative singularities are shown. The red solid line represents the source at a negative singularity, while the blue dashed line represents the source at a positive singularity. The experiment was conducted over 600 000 time-steps, with a total simulation time of 900 s, and other parameters identical to those in Fig. 4. The energy is calculated by summing the squared velocities of the simulated region. The magnitude of the energy in the figure has been appropriately scaled to highlight the relative variation of energy over time.
Two-layer model
In the second numerical experiment, we consider a two-layer medium model with a low-velocity upper layer. As shown in Fig. 11, the material properties of the two layers are as follows: the upper low-velocity layer has a P-wave velocity of 3.0 km s−1, an S-wave velocity of 1.2 km s−1 and a density of 1.8 g cm−3, while the lower high-velocity layer has a P-wave velocity of 6.0 km s−1, an S-wave velocity of 3.14 km s−1 and a density of 2.2 g cm−3. The thickness of the low-velocity layer is 3 km. The grid used in this model is similar to that employed in the first experiment, with one of the fine grid lines adjusted to align with the horizontal interface between the layers. This grid contains 501 076 grid points, representing a 44.79 per cent reduction compared to the 907 625 grid points required for a uniformly fine grid.
Snapshots of the vx and vy wavefields at t = 6.40 s in a two-layer medium model. The blue solid line indicates the position of the interface.
A single vertical force point source located at (x: 1793.28 m, y: −819.58 m) has source parameters of fc = 2.0 Hz and t0 = 0.7071 s. The smallest PPW occurs in the upper fine grid, where the grid spacing is 15.12 m, yielding a PPW of 15.76. This experiment also considers the influence of the free surface on seismic wave propagation. A free-surface condition is applied using a traction image method (Zhang & Chen 2006) at the top boundary, while PML boundary conditions are applied on the other three sides. The time-step is set to 0.001 s, with a total of 20 000 steps.
Nine receivers, represented by black triangles in Fig. 11, are placed at various locations throughout the model. Their detailed coordinates can be found in Table 1 under the ‘Experiment 2’ column. Receivers 1, 2 and 3 are located at the surface, while Receivers 6 and 9 are near the singularities. The dark blue solid line represents the boundary between the two layers.
As shown in Fig. 11, seismic waves reflect multiple times between the surface and the interface, with most energy concentrated in the low-velocity layer. This is evident in Fig. 12, where surface receivers record larger amplitudes than those in the lower high-velocity layer. The waves, generated by a point source (marked by a star), exhibit significant refraction and reflection at the layer boundary, creating complex wave patterns in the upper and lower layers. Additionally, Fig. 12 compares the results of CGFDM with GRTM, showing high consistency between our method and the reference solution. The quantitative errors, including EM and PM, are around 1–2 per cent, while EG and PG are close to 9.7–9.9, meeting the ‘excellent’ standard.
Comparison of the CGFDM and GRTM computed waveforms for the vx and vy components in the two-layer medium model.
Local low-velocity anomaly
In seismic studies, local low-velocity anomalies are sometimes encountered (Debayle et al. 2001; Han et al. 2021), requiring locally refined grids. To address this, we employed a grid refinement technique for a circular region (centred at x: 3780 m, y: −3780 m, with a radius of 1350 m) using four pairs of singularities arranged in a ring distribution, as shown in Fig. 13. This setup ensures that the grid is locally refined within the circular region, capturing the complex wave propagation behaviour with higher accuracy in this area. This grid consists of 344 305 grid points, representing a 70.21 per cent reduction compared to the 1 155 635 grid points required for a uniformly fine grid.
Wavefield snapshots of vx and vy at t = 1.40 s for the model containing a low-velocity anomaly. The dark blue line represents the boundary of the circular anomaly, and the green dots indicate the positions of the four pairs of singularities.
In the third numerical experiment, we use a point source (x: 1997.55 m, y: −1005.59 m) with fc = 8.0 Hz and t0 = 0.1768 s. The background medium has P-wave and S-wave velocities of 6.0 and 3.14 km s−1, respectively, and a density of 2.2 g cm−3. The medium properties are adjusted in the localized region, with P-wave and S-wave velocities of 3.0 and 1.2 km s−1, respectively, and a density of 1.8 g cm−3. The minimum PPW in this setup is 10.55. A time-step of 0.0005 s is used, with 8800 time-steps. Eight receivers are positioned at the locations marked by black triangles in Fig. 13, with Receivers 1, 2, 3 and 4 at the surface.
As illustrated in Fig. 13, the anomaly significantly impacts seismic wave propagation. Due to its much lower velocity than the surrounding medium, both P and S waves slow down and experience diffraction when passing through the anomaly, leading to notable wavefront distortion. Additionally, many scattering waves (as shown in Fig. 14) are generated due to multiple reflections within the anomaly. These scattering waves can be further amplified due to the interference and superposition of waveforms caused by reflections in this region.
Comparison of waveforms computed by CGFDM and FEM for the model containing a local low-velocity anomaly.
Due to the complexity of this model, we use the FEM as a reference solution. The nearly identical waveform comparison between the two methods is shown in Fig. 14, where the EM and PM fluctuate around 1 per cent to 2 per cent, while the EG values are above 9.6, and PG values are above 9.9. These quantitative comparisons indicate that the match between our method and the reference solution is at an ‘excellent’ level.
This simulation demonstrates the importance of locally refined grids when modelling wave propagation in the presence of subsurface heterogeneities. It provides a detailed approach to understanding how local anomalies influence seismic waves and how complex underground structures affect wavefields.
DISCUSSION AND CONCLUSION
In the first two experiments presented in this paper, the coarse-to-fine grid ratio is approximately 2:1, while in the third experiment, the ratio is 3:1. If larger ratios are needed, more pairs of positive and negative singularities can be added to the grid. The singularities used in this approach always appear in pairs, consisting of a positive and a negative singularity. While constructing difference templates at the singularity, we found that reusing one grid line at the negative singularity and removing one grid line at the positive singularity are arbitrary choices that do not affect the final results. Our method imposes no restrictions on the locations of the receivers, although there is a requirement for the source to be placed at least five grid points away from the singularities. Our method relies entirely on using difference templates and retains the computational efficiency of structured grids. The reduction in the number of grid points approximately reflects the reduction in computational cost. The third numerical experiment demonstrates the method's ability to accurately delineate subsurface anomaly boundaries during seismic wave simulations. This technique could potentially be applied to detect cavernous reservoirs and caves (Yang et al. 2012; Li et al., 2015) in future research.
This study successfully combines CGFDM with multiblock grid techniques for numerical simulations of seismic waves, achieving non-uniform grids by adding pairs of singularities. The treatment of singularities in our approach is extremely simple, maintaining the simplicity and efficiency of the FDM. The results from the three numerical experiments demonstrate a high level of agreement between our method and the reference solutions, and the method shows promising stability over extended simulation times, confirming the effectiveness of our approach in handling non-uniform grids. This study only considers the size of grid spacing, and future work will explore combining this approach with variable time-step methods. Although the examples presented are relatively simple, they are highly scalable, and the orthogonal boundaries of the grids allow for seamless integration, enabling the method to handle more complex seismic wave simulations. In three dimensions, implementing such non-uniform grids requires more sophisticated grid generation techniques, presenting certain challenges. We plan to explore these complexities in our future research.
DATA AVAILABILITY
The synthetic data are available from the corresponding author upon reasonable request.
The ICEM software can be downloaded at https://www.ansys.com/resource-center. The Pointwise software can be downloaded at https://www.cadence.com. The COMSOL Multiphysics software can be downloaded at https://www.comsol.com/acoustics-module.
ACKNOWLEDGEMENTS
We would like to express our sincere gratitude to Dr Peter Moczo, Dr Yder Masson and the editor, Dr Herve Chauris, for their invaluable feedback, insightful comments and diligent oversight during the review process. Their expertise and guidance have significantly enhanced the quality and clarity of the manuscript. This work is supported by the National Natural Science Foundation of China (42474065), Guangdong Provincial Key Laboratory of Geophysical High-resolution Imaging Technology (2022B1212010002) and High Level Special Funds (G03050K001). This work is supported by Center for Computational Science and Engineering at Southern University of Science and Technology. We thank Jialiang Wan for his help in programming the multigrid solver.
