Numerical simulation of P-wave propagation in 2D VTI media with a topographic free surface by a curvilinear-grid finite-difference method

Abstract

Topographic free surfaces and seismic anisotropy significantly modify seismic wave propagation. With increasing demands for high-precision subsurface imaging, advanced seismic prospecting and seismological imaging methods require simultaneous consideration of topographic free surface and anisotropy effects on P-wave propagation. However, conventional finite-difference methods struggle to accurately simulate both topographic free surface and anisotropy in first-order acoustic velocity-stress equations, failing to capture complex wavefield responses. To address this challenge, we extend the existing curvilinear-grid finite-difference method, originally developed for the elastic wave VTI system of velocity-stress equations, to the acoustic VTI velocity-stress system with topographic free surfaces. The primary difficulty lies in setting the shear-wave velocity to zero in the elastic wave equations, which causes the free-surface velocity-derivative constraint matrix for acoustic topographic free-surface models to become rank-deficient or severely ill-conditioned. Therefore, we derive the velocity-derivative constraint specific to the acoustic VTI system and use the Moore-Penrose generalized inverse to obtain the normal velocity derivatives needed for updating the free-surface boundary condition. The resulting formulation includes body-fitted curvilinear grids, collocated-grid DRP/opt MacCormack discretization, and traction-image boundary treatment. Numerical experiments on flat-surface, Gaussian-topography, and BP benchmark models demonstrate stable performance under the tested conditions, while near free-surface waveforms computed by our approach closely match independent reference solutions. This method provides a reliable forward-modeling foundation for 2D acoustic-wave applications that require simultaneous, high-accuracy consideration of irregular topography and anisotropy.

1 Introduction

Numerical simulation of seismic wave propagation in complex media is a foundational technology for investigating wave propagation laws, unraveling complex wavefield mechanisms, and enabling high-precision seismic imaging and inversion. With the rapid advancement of wave-equation-based methods such as full-waveform inversion (e.g., Tarantola, 1984; Virieux and Operto, 2009) and reverse-time migration (e.g., Baysal et al., 1983), seismic wave simulations are required to account for the comprehensive effects of the topographic free surface and medium anisotropy on wavefields. The finite difference method (FDM; e.g., Virieux, 1984; 1986; Levander, 1988; Graves, 1996; Igel et al., 1995; Saenger and Bohlen, 2004; Zhu et al., 2009; Lisitsa and Vishnevskiy, 2010), finite-element method (e.g., Zhu and Dorman, 2000), spectral-element method (SEM; e.g., Komatitsch and Tromp, 1999; Komatitsch et al., 2000), and boundary-element method (e.g., Brebbia, 1981; Dominguez, 1993; Ge et al., 2005) have been widely applied in seismic wavefield simulations. Among these methods, the FDM is widely applied in seismological and exploration geophysics research due to its algorithmic simplicity, high computational efficiency, ease of parallel implementation and lower computing cost, establishing it as a widely adopted numerical simulation technique in industrial geophysical applications.

The free surface is a strong discontinuity in seismic wave propagation, governing numerous wave propagation phenomena. The topographic free surface directly affects wave-propagation behavior and induces multiscale wavefield effects such as scattering, diffraction, focusing, and defocusing. These effects distort the originally clear reflected wavefields and systematically alter the amplitude, phase, and travel-time characteristics of seismic records. For seismic imaging and inversion, failure to accurately characterize topographic free surfaces numerically may lead to errors in wave propagation paths, false structural imaging, and distortions in the local amplitude estimates. Therefore, in high-precision seismic wave numerical simulations, topographic free surface effects must be treated as independent boundary conditions with refined processing, rather than simply approximated as minor perturbations of a horizontal free surface. The influence of an undulating free surface is particularly important for seismic prospecting and seismological imaging methods. Rugged topography changes source and receiver elevations, local incidence and reflection angles, and near-surface propagation paths; it also generates scattering, diffraction, and spatially varying focusing and defocusing. If these effects are omitted from the forward operator, the resulting phase and amplitude mismatches may cause cycle skipping, false reflectors, and biased model updates. Previous studies have demonstrated that irregular topography must be accurately represented in elastic full-waveform inversion and reverse-time migration to avoid waveform and imaging errors (e.g., Borisov et al., 2018; Liu et al., 2017). Attenuation and VTI anisotropy are also important physical factors that should be consistently incorporated into wave-equation-based inversion (e.g., Song et al., 2023; Chen et al., 2026).

Conventional regular Cartesian grids in finite-difference numerical simulations struggle to accurately conform to undulating free surfaces. While the commonly used staircase approximation methods are simple to implement, they discretize the topographic free surface into stepped boundaries, introducing non-physical scattering, numerical dispersion, and phase errors. To address this issue, researchers have developed complex-terrain handling strategies, such as the air layer or vacuum layer method (e.g., Boore, 1972; Graves, 1996), the boundary-corrected difference method (e.g., Jih et al., 1988; Bayliss et al., 1986; Kristek et al., 2002), the immersed boundary method (e.g., Li et al., 2020), and the curvilinear grid method (e.g., Fornberg, 1988; Hestholm and Ruud, 1994; 1998; Zhang and Chen, 2006; Wang and Liu, 2007; Appelö and Petersson, 2009; Lan and Zhang, 2011; Zhang, Zhang, and Chen, 2012; Zhang et al., 2014). Among these methods, the curvilinear grid method maps an undulating ground surface in physical space to regular boundaries in computational space through coordinate transformation. This approach maintains the efficiency of structured grids while better conforming to complex topographies, avoids strong scattering issues caused by staircase approximation, and facilitates the relatively easy application of free-surface boundary conditions (e.g., Hestholm and Ruud, 1998; Zhang and Chen, 2006; Zhang et al., 2012; Sun et al., 2016).

Besides the topographic free surface, seismic anisotropy represents another critical factor influencing seismic wave propagation. Actual crustal media, particularly sedimentary rock formations, shale, thin interbedded layers, and rocks with oriented fractures, generally do not satisfy the isotropic assumption (e.g., Alkhalifah, 1998). Thomsen (1986) noted that most rock media exhibit weak to moderate elastic anisotropy. Even weak anisotropy may exert non-negligible impacts on seismic wave propagation, reflection travel-time relationships, and velocity analysis in exploration seismology. Transversely isotropic media are commonly used anisotropic models in seismic exploration. Among these, VTI media can be characterized by vertical P- and S-wave velocities and Thomsen weak anisotropic parameters. Therefore, it is essential and feasible to simultaneously account for both topographic free surfaces and seismic anisotropy in high-precision numerical simulations. In seismic exploration, the primary research objective is typically to obtain the subsurface P-wave velocity structure, which supports reservoir-structural interpretation, lithology identification, and subsequent imaging inversion. Compared with full elastic-wave simulation, acoustic or pseudo-acoustic equations can suppress shear wave and surface wave components while preserving the primary P-wave propagation characteristics, thereby reducing the computational cost (e.g., Alkhalifah, 2014; Xu et al., 2020). In actual surface observations, shear waves and surface waves are often regarded as interference components, and shear wave energy attenuates rapidly in near-surface complex media. Therefore, under the premise of balancing computational efficiency and P-wave simulation accuracy, anisotropic acoustic equations are widely used in applications such as migration imaging, reverse-time migration, and full waveform inversion. However, most existing studies on anisotropic acoustic simulation are based on horizontal free surfaces or regular grids, and insufficient attention has been paid to P-wave propagation under the combined action of undulating free surfaces and medium anisotropy.

Commonly used finite difference schemes include staggered-grid finite difference, rotated staggered-grid finite difference, and collocated-grid finite difference. The conventional staggered-grid finite difference method leverages the symmetric structure of the isotropic elastic wave equation to efficiently simulate seismic wave propagation under horizontal surface conditions. However, when topographic relief is significant, traditional staggered grids become complex in handling mixed derivatives under curvilinear coordinates, boundary conditions, and variable position matching. On the other hand, collocated-grid methods feature simple variable arrangements, facilitating unified treatment of spatial derivatives and boundary conditions in curvilinear coordinates. Nevertheless, conventional central differences on collocated grids tend to produce odd-even decoupling and unphysical high-frequency oscillations. Therefore, it is necessary to adopt high-precision difference schemes with low-dispersion and low-dissipation characteristics to improve computational stability and wavefield simulation accuracy. The DRP/opt MacCormack scheme (e.g., MacCormack, 1969; Hixon, 1997; Hu et al., 1996) alternately uses forward and backward one-sided difference operators. While maintaining high-order accuracy, it implicitly introduces appropriate numerical dissipation, which helps suppress unphysical oscillations in collocated grids. Thus, it is suitable for solving first-order velocity-stress wave equations under curvilinear grids (e.g., Zhang et al., 2012; Sun et al., 2016). Although Sun et al. (2016) applied the curvilinear grid-based traction imaging method to 2D anisotropic elastic waves, they did not extend its application to anisotropic acoustic wave equations. While deriving an acoustic wave equation by setting the shear wave velocity to zero in elastic wave formulations may seem conceptually straightforward, none of these authors attempted or proposed such an extension. Our research identifies a critical reason for this oversight: when the shear wave velocity is set to zero, the boundary condition implementation in their elastic wave methods becomes numerically unstable near the free surface for acoustic wave equations, failing to produce reliable solutions. This numerical instability constitutes the key problem addressed in the present study. Furthermore, simultaneously accounting for terrain undulation and VTI anisotropy is not a simple combination of two separate effects. The coordinate transformation introduces distinct metric coefficients that couple derivatives along computational directions, while the VTI constitutive equation links horizontal and vertical stress responses. Consequently, the normal-velocity-derivative constraint matrix depends on both the local surface geometry and anisotropic parameters, and may become rank-deficient or severely ill-conditioned. Therefore, the elastic traction imaging method cannot be directly applied to the pseudo-acoustic VTI system without specialized boundary condition modifications.

In this research, we extend the existing curvilinear-grid finite-difference model for elastic wave equations to acoustic wave equations in 2D anisotropic VTI media with a topographic free surface. The main contribution is the formulation of a velocity-derivative constraint and a Moore-Penrose generalized inverse closure for the normal derivatives needed for stably implementing the free surface boundary conditions in VTI media. This approach enables the application of the traction image method under the condition that the boundary matrix is rank deficient or severely ill-conditioned. Our numerical experiments comprehensively examine the combined effects of topography, VTI anisotropy, and shallow low-velocity layers, incorporating both free-surface boundary conditions and absorbing layers. The calculated wavefield snapshots, common-shot gathers and single trace waveforms are employed to evaluate both the numerical stability and the physical processes affecting P-wave travel times, amplitudes, scattering and energy redistribution in the near surface region.

2 Methodology

2.1 The Hooke’s law for acoustic wave in 2D VTI media

In this study, we adopt the widely used Hooke’s law for VTI anisotropic media, as proposed in the classic work of Alkhalifah (1998). Using the Voigt notation, Hooke’s law for elastic waves in VTI media can be expressed aswhere represents components of the stress tensor, denotes components of the strain tensor, and indicates components of the elastic stiffness tensor. The relationship between the five independent parameters of VTI media (, , , , ) and the vertical P-wave velocity (), S-wave velocity (), and Thomsen parameters (, , ; Thomsen, 1986) can be expressed as , , , , and . Hooke’s law for acoustic waves in VTI media can be derived by setting the shear-wave velocity to zero () in Equation 1. This approximation, discussed by Alkhalifah (1999), Alkhalifah (2000) and Grechka et al. (2004), leads to the following form:

Hooke’s law for acoustic VTI media, presented in Equation 2, can be reduced to two stress components (Duveneck and Bakker, 2011). Here, the horizontal stress component is defined as , and the vertical stress component is denoted as . Accordingly, the first-order velocity-stress formulation of the acoustic wave equation for VTI media can be expressed in 2D Cartesian coordinates by assuming spatial derivatives with respect to the y-axis are zero:

This set of equations (Equations 3–6) is widely employed in studies of anisotropic acoustic wave propagation. Since the stress component is absent, the free-surface boundary condition for a planar surface with normal vector along the z-axis simplifies to

2.2 The velocity-stress equation in curvilinear coordinates

For the topographic free surface, we discretize the model using curvilinear grids to generate body-fitted grids conforming to the topographic free surface. Through the coordinate transformation shown in Figure 1, body-fitted grids map curved surfaces in physical space (-) to uniform grids in computational space (-), making them suitable for finite-difference numerical calculations. Detailed discussions on the coordinate mapping method can be found in Zhang and Chen (2006).

FIGURE 1

The acoustic wave equations in Equations 3–6 can be transformed into curvilinear coordinates to accommodate body-fitted grids. In curvilinear coordinates, the acoustic VTI wave equation can be expressed using the chain rule in Equations 7, 8:and expressed in the following form in Equations 9–12:where the coordinate mapping coefficients , , , and can be directly derived from the Cartesian coordinates of the curvilinear grid points. The detailed method for these calculations is described in Zhang and Chen (2006).

2.3 The DRP/opt MacCormack finite difference method

The collocated-grid central difference scheme exhibits numerical instability known as odd-even decoupling when solving wave equations. To suppress this phenomenon, dissipation or filtering terms must be incorporated into the difference scheme (e.g., Zhang and Chen, 2006). The MacCormack difference scheme implicitly introduces high-frequency dissipation through forward and backward one-sided differences (e.g., MacCormack, 1969; Zhang and Chen, 2006; Sun et al., 2016). For convenience of discussion, consider the following 1D elastic wave equation:where denotes the coefficient matrix and represents the wavefield vector. We adopt the DRP/opt MacCormack difference scheme proposed by Hixon (1997) to achieve fourth-order spatial accuracy and high grid resolution. In the one-dimensional case, it can be expressed in Equations 14, 15 as:where, and represent the forward and backward 5-point stencil one-sided differences at grid point , respectively, with difference coefficients , , , and . Numerical experiments demonstrate that the 4th-order Runge-Kutta scheme is sufficient to ensure the required accuracy for seismic wave propagation simulations (Zhang et al., 2012; Sun et al., 2016). Therefore, the fourth-order Runge-Kutta scheme for the 1D elastic equation (Equation 13) can be expressed in Equations 16–20 as:where denotes the time step, , , and the coefficients , , , , , , , . For discussion convenience, the above Runge-Kutta scheme can be simplified in Equation 21 as:

Following the 1D notation, the right-hand sides of Equations 3–6 are computed using the combined MacCormack operators below in Equations 22–25, where subscripts and denote computational-coordinate directions:

Following Sun et al. (2016), the fourth-order Runge-Kutta scheme in Equations 26–29 alternates the MacCormack directional combinations over four consecutive time steps to reduce directional bias:

While minor discrepancies exist between Equations 26–29 and the formulations presented by Sun et al. (2016), our numerical tests demonstrate that this time-evolution scheme exhibits robust stability. For comprehensive details regarding this approach, readers are referred to Zhang and Shen (2010) and Sun et al. (2016).

The stability condition of the 2D DRP/opt MacCormack scheme is expressed in Equation 30 as:where represents the minimum spatial grid spacing, and denotes the maximum velocity in the medium.

2.4 The free-surface boundary condition

To ensure generality, the free-surface boundary condition is imposed in the 2D computational space, with the assumption that the -axis of the computational space always coincides with the free-surface topography. In this case, the normal vector of the free-surface is given by , where and are unit vectors along the and axes in the Cartesian coordinate system, and . According to the free-surface boundary condition, the traction force must satisfy , leading to the following expression in Equations 31, 32:

Similar to free-surface boundary condition handling methods for the elastic wave equation (e.g., Zhang and Chen, 2006; Sun et al., 2016), the acoustic wave equation can also be used to derive the constraint form of the free-surface boundary condition on velocity derivatives based on Hooke’s law as follows:where the coefficient matrices are:

Equations 33–35 are formal solutions derived by analogy with the elastic wave equation. The matrix requires inversion and its determinant is . However, in practical applications, may become singular or severely ill-conditioned under specific combinations of topographic surface geometry and anisotropic parameters, especially at locally horizontal or vertical surface segments and in the elliptically anisotropic limit , including the isotropic case. This rank deficiency is associated with the assumption of zero shear-wave velocity in acoustic wave equations, which leads to the disappearance of the component and reduces the completeness of the free surface boundary constraints. Through extensive numerical tests, we ultimately adopted the Moore-Penrose generalized inverse matrix method (e.g., Moore, 1920; Penrose, 1955) to replace with in the formal full-rank solution, and found that this method ensures stable algorithmic operation. Mathematically, represents the minimum-norm solution of the equation, which guarantees continuity and stability of the computation process.

3 Numerical experiments

To validate the stability of the proposed algorithm and its capability to handle topographic free surfaces, seismic anisotropy, and shallow low-velocity layers, we employ three distinct models—a horizontally layered model, a topographic free-surface model, and the widely used BP model in the oil industry—to conduct numerical experiments using the proposed curvilinear-grid finite-difference method. In these experiments, we investigate the effects of topographic free-surfaces, seismic anisotropy, and low-velocity layers on wavefield snapshots, shot gathers, and theoretical waveforms.

3.1 The flat free-surface model

The horizontal free-surface model established in this study spans horizontally from x = 37.5 km to x = 49.5 km and vertically from z = −4.5 km to z = 0.0 km, where the free surface lies at z = 0.0 km. The P-wave velocity models are illustrated in Figures 2A,B. The model for Figure 2B is obtained from Figure 2A by adding a near-surface low-velocity layer. The anisotropic parameter models for δ and ε are shown in Figures 2C,D. Two sets of models were designed for this study. The first set uses the velocity model without the low-velocity layer (Figure 2A) to calculate acoustic wavefields in both isotropic and anisotropic media. By contrast, the second set uses the velocity model containing the low-velocity layer (Figure 2B) to calculate acoustic wavefields in isotropic and anisotropic media. Both sets of models adopt a grid spacing of 5.0 m and a time-step of 0.1 ms. A total of 900 receivers are deployed along the free-surface x-axis, starting at x = 39.010 km with a 10.0-m spacing between adjacent receivers. The explosive source is placed at x = 43.505 km, 10.0 m beneath the free surface (2 grid points away from the surface), and a 40.0-Hz Ricker wavelet serves as its source-time function. Figures 3, 4 present the numerical simulation outcomes from the first set, and Figures 5, 6 present those from the second set.

FIGURE 2

FIGURE 3

FIGURE 4

FIGURE 5

FIGURE 6

In Figure 3, we present snapshots of the Vz-component wavefields at 0.8 s, 1.0 s and 1.2 s, which are calculated from the isotropic and VTI models built upon the P-wave velocity model in Figure 2A. A comparison between results for the isotropic media model (Figures 3A,C,E) and the VTI media model (Figures 3B,D,F) clearly reveals reflected waves generated at the free surface and internal interfaces. Wavefields propagate properly around the free surface and internal interfaces, indicating that our algorithm maintains stable wavefield computation under the simplified condition of a horizontal free-surface and isotropic media. Three blue vertical dashed lines are placed at x = 39.5 km, 40.5 km and 41.5 km for comparing the horizontal propagation characteristics of wavefields. Comparison between Figures 3C,D shows that the wavefront in the VTI media has propagated past x = 40.5 km, while the wavefront in Figure 3C has just reached x = 40.5 km. A similar phenomenon can be observed in Figures 3E,F. This indicates that the wavefields in the VTI media propagate faster along the horizontal direction, which is consistent with the δ and ε parameter models shown in Figures 2C,D. To further investigate wavefields near the free surface, we present in Figure 4 the common-shot gathers recorded by free-surface receivers and waveform comparisons of the trace at an offset of −2.0 km. From the common-shot gathers in Figures 4A,B, we observe that the VTI media model exhibits similar reflection phases, but the direct wave in the isotropic model is stronger than that in the VTI model, while the 2^nd-layer reflected wave in the VTI model has higher amplitude. At far offsets, the travel time of waves in VTI is shorter than that in isotropic media. The same conclusion can be drawn from the waveform comparison at an offset of −2.0 km.

In Figure 5, we present snapshots of the Vz-component wavefields at 0.8 s, 1.0 s and 1.2 s, calculated using both isotropic and VTI models with a low-velocity layer shown in Figure 2B. The wavefield snapshots in Figure 5 are significantly more complex than those in Figure 3, which is attributed to the presence of the shallow low-velocity layer. By comparing the isotropic media results (Figures 5A,C,E) with the VTI media results (Figures 5B,D,F), we can clearly observe the complex reflected waves from the free surface and internal interfaces. The wavefields propagate normally near the free surface and internal interfaces, indicating that our algorithm maintains stable performance even in the presence of a shallow low-velocity layer. By comparing Figures 5C,D, we observe that the wavefront in VTI media has propagated past x = 40.5 km, while the wavefront in the isotropic model has not yet reached x = 40.5 km. A similar phenomenon is observed in Figures 5E,F. This indicates that wavefields in VTI media propagate faster along the horizontal direction, consistent with the δ and ε parameter models shown in Figures 2C,D. To further investigate wavefields near the free surface, we present in Figure 6 the common-shot gathers recorded by surface receivers and waveform comparison of a trace at an offset of −2.0 km. From the common-shot gathers in Figures 6A,B, we observe that the VTI media model exhibits similar reflection phases to the isotropic media, whereas the direct wave in the isotropic model shown in Figure 6C is slightly stronger than that in the VTI model, and the reflected wave amplitudes are comparable between the two models. At far offsets, the travel time of waves in VTI is shorter than that in isotropic media.

From the horizontal model experiments, we observe that anisotropy affects the travel times and amplitudes of direct and reflected waves, with the travel time effect becoming more pronounced as the offset increases. The shallow low-velocity layer reduces the amplitude differences of direct and reflected waves between isotropic and anisotropic media, but does not affect the travel time differences at far offsets. These horizontal model tests also demonstrate that our algorithm can stably and accurately simulate seismic wave propagation in anisotropic media.

3.2 Topographic free-surface model

The topographic free-surface model established in this study is modified from the horizontal free-surface model shown in Figure 2. Specifically, the horizontal free surface in Figure 2 is replaced with a Gaussian hill, whose mathematical expression is given by: . The protruding part of the Gaussian hill adopts the same medium parameters as the surrounding single layer. The P-wave velocity models are illustrated in Figures 7A,B, where the model in Figure 7B is derived from Figure 7A by adding a near-surface low-velocity layer. The anisotropic parameter models for δ and ε are presented in Figures 7C,D. Two sets of numerical simulation experiments were designed in this study. The first set uses the velocity model without the low-velocity layer (Figure 7A) to simulate acoustic wavefields in both isotropic and anisotropic media. The second set uses the velocity model with the near-surface low-velocity layer (Figure 7B) to simulate acoustic wavefields in both isotropic and anisotropic media. Both model configurations adopt a grid spacing of 5.0 m and a time interval of 0.1 ms. A total of 900 receivers are deployed along the x-axis at the free surface, starting from x = 39.010 km with an interval of 10.0 m between adjacent receivers. The explosive source is placed at x = 43.505 km, 10.0 m underground (two grid points away from the surface), with a 40.0 Hz Ricker wavelet as the source time function. Figures 8, 9 present the results of numerical simulation experiments using the first model configuration, while Figures 10, 11 show the results of experiments using the second model configuration.

FIGURE 7

FIGURE 8

FIGURE 9

FIGURE 10

FIGURE 11

In Figure 8, we present snapshots of the Vz-component wavefields at 0.8 s, 1.0 s, and 1.2 s, calculated using both isotropic and VTI models with the Vp model shown in Figure 7A. By comparing the isotropic media results (Figures 8A,C,E) with the VTI media results (Figures 8B,D,F), we can clearly observe reflected waves from the free surface and internal interfaces. The wavefields propagate normally near the free surface and internal interfaces, indicating that our algorithm can stably compute wavefields even when the model degenerates to an isotropic model with a topographic free surface. By comparing Figures 3C,D, we observe that the wavefront in the VTI medium is closer to x = 40.5 km than that in the isotropic medium. Since the only difference between the models in Figure 3 and the current model is the topography, this indicates that topographic relief affects wave travel times. To further investigate wavefields near the free surface, we present in Figure 9 the common-shot gathers recorded by surface receivers and a waveform comparison of a trace at an offset of −2.0 km. From the common-shot gathers in Figures 9A,B, we observe that the VTI media model exhibits similar reflection phases to the isotropic model. In Figure 9C, the direct wave in the isotropic model matches well with that in the VTI model, while the second-layer reflected wave in the VTI model has a higher amplitude. The first layer in the VTI model shows weak anisotropy, which is consistent with field observations. This is because topographic relief causes thickening of the isotropic layer, so the direct wave is not affected by the anisotropic layer, resulting in excellent waveform matching with the isotropic model. At far offsets, the travel time of waves in the VTI medium is shorter than that in the isotropic medium.

In Figure 10, we present snapshots of the Vz-component wavefields at 0.8 s, 1.0 s, and 1.2 s, calculated using both isotropic and VTI models that include the low-velocity layer shown in Figure 7B. The wavefield snapshots in Figure 10 are significantly more complex than those in Figure 8, which can be attributed to the presence of the shallow low-velocity layer. By comparing the isotropic media results (Figures 10A,C,E) with the VTI media results (Figures 10B,D,F), we can clearly observe complex reflected waves generated from the free surface and internal interfaces. The wavefields propagate normally near the free surface and internal interfaces, indicating that our algorithm maintains stable performance even when simulating wave propagation in the presence of a shallow low-velocity layer. To further investigate wavefields near the free surface, we present in Figure 11 the common-shot gathers recorded by surface receivers and a waveform comparison of a trace at an offset of −2.0 km. From the common-shot gathers in Figures 11A,B, we observe that at far offsets, the amplitudes of direct waves are weaker than those of reflected waves in both isotropic and VTI media. This is significantly different from the common-shot gathers in Figure 9, which is caused by the shallow low-velocity layer concentrating the energy excited by the source in the near-offset region. From the waveform comparison in Figure 11C, due to the strong anisotropy in the deep part of the model, the travel time difference of deep reflected waves between the VTI medium and the isotropic medium is significant.

From the topographic free-surface model experiments, we observe that seismic anisotropy and free-surface topography affect the travel times and amplitudes of direct and reflected waves. The combination of the shallow low-velocity layer and topography concentrates seismic energy in specific regions, leading to variations in the amplitudes of direct and reflected waves, but does not affect the travel time differences in the deep part of the model. These topographic free-surface model experiments also demonstrate that our algorithm can stably and accurately simulate seismic wave propagation in anisotropic media.

3.3 The absorbing layer and free surface boundary conditions

In seismic exploration and imaging studies, the free-surface boundary is sometimes replaced with absorbing boundary layers during wavefield analysis and shot gather processing to eliminate surface reflections. Our proposed method supports both the perfectly matched layer (PML) technique (Zhang and Shen, 2010) and the nonreflecting layer approach introduced by Cerjan et al. (1985). The Cerjan-type absorbing boundary is adopted in the current implementation due to its straightforward implementation and numerical stability. To demonstrate this capability, we removed the free-surface boundary condition from both the flat-topography models (Figures 5B, 6B) and the topographic models (Figures 10B, 11B), and extended the computational domain upward by 24 grid points to incorporate absorbing layer boundary conditions. New wavefield snapshots (Figures 12A,C) and shot gathers (Figures 13A,C) were generated using this modified setup. For comparative analysis, the corresponding wavefield snapshots (Figures 12B,D) and shot gathers (Figures 13B,D) were obtained from the original free-surface models presented in Figures 5B, 6B, 10B and 11B.

FIGURE 12

FIGURE 13

Figure 12 displays a 2 × 2 comparison of seismic wavefield snapshots showing the vertical velocity component (Vz) in the VTI media at 0.80 s, with all subplots sharing identical axes and a unified color bar. Figure 12A presents the wavefield under flat topography with a Cerjan-type absorbing boundary, where seismic energy propagates in a clear hemispherical pattern concentrated within 0.0–3.0 km depth, with no visible surface reflections as waves reaching the top are fully attenuated. In contrast, Figure 12B features a free surface boundary on flat terrain, producing distinct reflected wave trains that form concentric ring structures and create complex interference patterns with the incident wavefield, leading to disordered energy distribution near the surface. For undulating topography models, Figure 12C demonstrates that the absorbing boundary effectively suppresses terrain-induced interference, maintaining a relatively undisturbed hemispherical wavefield despite the irregular surface, with no upward-propagating reflected waves observed. Meanwhile, Figure 12D shows that a free surface boundary on undulating terrain generates strong, irregularly propagating reflected waves, which interfere severely with the incident wavefield and cause significant distortion of the wavefield morphology beneath the surface. These comparisons clearly validate the effectiveness of absorbing boundaries in eliminating surface reflections and preserving wavefield integrity, particularly in complex topographic conditions, while highlighting the substantial wavefield interference caused by free surface reflections. These conclusions are further corroborated by their corresponding shot gathers shown in Figure 13.

3.4 Synthetic waveform comparisons with SEM solutions

To quantitatively validate the accuracy of the proposed FDM for solving 2D acoustic velocity-stress equations with flat and topographic free surfaces, we use waveforms computed by the SEM package SPECFEM2D (https://github.com/SPECFEM/specfem2d) as reference solutions. We constructed two acoustic models: one with a flat topography and one with a Gaussian hill topography. In each model, five vertically distributed receivers were deployed from the free surface to the medium interior (Figure 14). Receivers R1-R5 correspond to the flat free-surface model, while receivers R6-R10 are associated with the topographic free-surface model. The recommended accuracy for SPECFEM2D is 5 elements per wavelength. To ensure high-precision waveforms, we increased this parameter to 15 elements per wavelength. This configuration produced reference waveforms at the five receiver positions for both the flat model (Figure 15) and the undulating model (Figure 16). As shown in Figures 15, 16, the Vz and Vx waveforms computed by the proposed FDM exhibit excellent agreement with the high-precision SEM waveforms, with cross-correlation coefficients exceeding 0.99. The number of grid points per wavelength for the proposed FDM computing is approximately 8.0. These results confirm that our method can accurately compute wavefields in both flat and undulating topographic models.

FIGURE 14

FIGURE 15

FIGURE 16

3.5 The BP model

In this experiment, we use the industry-standard BP model (Gray and Marfurt, 1995) for testing to demonstrate two key capabilities of our algorithm: stable wavefield computation in complex isotropic and anisotropic models, and the ability to quantify the combined effects of anisotropy and free-surface topography on acoustic wavefields. Due to the large scale of the full BP model, we selected a representative subregion for our computational experiments, spanning horizontally from x = 37.5 km to x = 49.5 km and vertically from z = −2.0 km to approximately z = 4.0 km at the free surface. The P-wave velocity model is illustrated in Figure 17A. The anisotropic parameters δ and ε, along with the density model, were derived from the P-wave velocity. The anisotropic-parameter models for δ and ε are shown in Figures 17C,D, respectively, while the density model is illustrated in Figure 17B. We performed wavefield computations for both isotropic and VTI models using these parameters, with a consistent grid spacing of 5.0 m and a time interval of 0.1 ms implemented in both simulations. A total of 900 receivers are deployed along the topographic free surface starting from x = 39.010 km, with an interval of 10.0 m between adjacent receivers. An explosive source was placed at x = 43.625 km, 15.0 m below the free surface, using a 40.0 Hz Ricker wavelet as the source time function. The numerical simulation results are presented in Figures 18, 19, which demonstrate the robustness of our algorithm in handling complex free-surface geometry and anisotropic properties.

FIGURE 17

FIGURE 18

FIGURE 19

In Figure 18, we present snapshots of the Vz-component wavefields at 0.8 s, 1.0 s, and 1.2 s, computed using both isotropic and VTI models based on the P-wave velocity model shown in Figure 17A. By comparing results from isotropic media (Figures 18A,C,E) with those from VTI media (Figures 18B,D,F), we clearly observe reflected waves from the free surface and internal interfaces. The wavefields propagate uniformly near the free surface and internal interfaces, demonstrating that our algorithm can stably compute wavefields in complex media models with topographic free surfaces. By comparing Figures 18C,D, we observe that the wavefront in the VTI media has propagated beyond x = 40.5 km, while the wavefront in Figure 18C has just reached x = 40.5 km. This difference indicates that seismic anisotropy affects wave travel times. To further investigate wavefield behavior near the free surface, we present in Figure 19 the common-shot gathers recorded by surface receivers and a waveform comparison of a single trace at an offset of −2.0 km. From the common-shot gathers in Figures 19A,B, we observe that the VTI media model exhibits similar reflection phases to the isotropic model. In Figure 19C, the direct wave in the isotropic model shows distinct amplitude and phase differences compared to that in the VTI model. Additionally, the VTI model displays some high-amplitude reflection phases, which can be attributed to seismic anisotropy. At far offsets, wave travel times in the VTI medium are shorter than those in the isotropic medium.

For the BP model experiment, we find that the effects of seismic anisotropy on the amplitude and phase of direct and reflected waves cannot be ignored, and topographic variations cause changes in wavefield propagation. Therefore, these factors need to be considered in seismic imaging processes. This numerical simulation experiment also demonstrates that our proposed algorithm can stably and accurately compute acoustic wave propagation in complex anisotropic media with topographic free surfaces.

4 Discussion

The numerical experiments presented in this study demonstrate the effectiveness of the proposed curvilinear-grid finite-difference method for simulating P-wave propagation in 2D VTI media with topographic free surfaces. The key innovation, using the Moore-Penrose generalized inverse to address rank deficiency in the free-surface velocity-derivative constraint matrix, enables stable and accurate wavefield modeling where conventional elastic-based traction image methods fail. This section discusses the implications of our findings, limitations of the current approach, and directions for future research.

4.1 Physical implications of anisotropy and topography

Our results consistently show that VTI anisotropy accelerates horizontal wave propagation relative to isotropic media, as evidenced by earlier wavefront arrival at far offsets in both flat and topographic models (Figures 3,5,8,10,18, , , , ). This aligns with theoretical predictions from Thomsen parameters (Thomsen, 1986), where positive δ and ε values enhance horizontal velocity. The amplitude differences between isotropic and VTI wavefields (e.g., Figures 4C, 6C, 9C) highlight the anisotropic influence on energy partitioning, which could impact amplitude-based inversion workflows if neglected.

Topographic free surfaces introduce complex wavefield phenomena, including scattering, diffraction, and focusing (Figures 8, 10, 12D). When combined with shallow low-velocity layers, topography concentrates seismic energy in near-offset regions (Figure 11), altering amplitude distributions and potentially biasing velocity estimates in conventional processing. The BP model experiments further demonstrate that topographic variations and anisotropy together distort wavefronts and reflection phases (Figure 18), emphasizing the need for their joint consideration in high-precision imaging.

4.2 Numerical stability and accuracy

The stability of the proposed method is verified across diverse scenarios: flat surfaces, Gaussian topography, and the complex BP model. The DRP/opt MacCormack scheme effectively suppresses collocated-grid oscillations, while the Moore-Penrose generalized inverse ensures boundary condition stability even for rank-deficient matrices (e.g., at horizontal surface segments or ε = δ). Waveform comparisons with SEM solutions (Figures 15, 16) confirm high accuracy, with cross-correlation coefficients exceeding 0.99, validating the method’s potential for quantitative applications like full-waveform inversion. Notably, the Cerjan-type absorbing boundary effectively eliminates surface reflections (Figures 12, 13), providing flexibility for applications where free-surface multiples are considered noise.

4.3 Limitations and future work

Despite its demonstrated strengths, the proposed method has inherent limitations that require consideration. First, the current implementation is restricted to 2D models. Extending the framework to 3D scenarios would necessitate addressing more complex curvilinear coordinate transformations and formulating 3D VTI constitutive relations, both of which introduce significant numerical challenges related to computational stability, memory management, and algorithm scalability. Second, curvilinear grids combined with high-order finite-difference schemes inherently increase computational demands relative to conventional Cartesian-grid methods. While optimizations such as adaptive mesh refinement or GPU acceleration could mitigate this overhead by concentrating computational resources on critical wavefield regions, further efficiency improvements remain a key area for development.

We acknowledge that computational efficiency is highly contingent on program optimization strategies and rigorous implementation practices, which can vary across different numerical methods. To establish a comprehensive understanding of our method’s performance relative to existing techniques, we plan to benchmark against state-of-the-art wave equation migration algorithms. For example, the one-step sparse-matrix method for full wave equation depth migration proposed by Li et al. (2024) demonstrates that sparse matrix multiplication can reduce computational costs while maintaining imaging accuracy, which offers valuable insights for optimizing our high-order finite-difference framework for full wave equation depth migration. Additionally, the generalized two-way phase shift operator for anisotropic media developed by Li et al. (2025) provides a promising approach for handling VTI wavefield propagation without high-order corrections, which could inform our future efforts to extend the method to complex anisotropic environments.

To ensure a fair and robust comparison, we plan to collaborate with experts specializing in advanced migration implementations to conduct a systematic quantitative analysis in future work. This analysis will account for factors such as parallelization strategies, memory usage, and problem-specific scaling behavior, while also exploring synergies with complementary techniques like sparse matrix acceleration and generalized phase shift operators.

5 Conclusion

In this study, we extend the curvilinear-grid finite-difference method to simulate acoustic wavefield propagation in anisotropic media with topographic free surfaces. The traction image method, originally developed for elastic wave equations, cannot be directly applied to acoustic wave equations due to the absence of a shear modulus. This absence may render the coefficient matrix governing the free-surface boundary condition singular or severely ill-conditioned for specific surface geometries and anisotropic parameters, making direct inversion inapplicable at those points. Our proposed algorithm addresses this issue by employing the Moore-Penrose generalized inverse method to obtain the minimum-norm closure of the rank-deficient boundary system. We validate the algorithm through numerical tests on three representative models: a planar free surface model, a topographic free surface model, and the complex BP model with topographic free surfaces. Experimental results demonstrate that our method can stably compute acoustic wavefields in all three models, confirming its robustness across different geological scenarios. While near free-surface waveforms computed by our approach closely match independent reference solutions. This method delivers a reliable forward-modeling foundation for 2D VTI acoustic wave applications that require simultaneous, high-accuracy consideration of irregular topography and anisotropy. The principal improvement of this work is that the Moore-Penrose closure enables the traction image free-surface treatment to be extended from elastic equations to the acoustic VTI system, while retaining the body-fitted curvilinear grid and the low-dispersion collocated-grid formulation.

References

• 1

  AlkhalifahT. (1998). Acoustic approximations for processing in transversely isotropic media. Geophysics63 (2), 623–631. 10.1190/1.1444361

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AlkhalifahT. (1999). An acoustic wave equation for Anisotropic Media. 69th Ann. Intern. Mtg., Soc. Expl. Geophys., Expand. Abstr., 722–725.

  • Google Scholar
  • View reference in article

• 3

  AlkhalifahT. (2000). An acoustic wave equation for anisotropic media. Geophysics65 (4), 1239–1250. 10.1190/1.1444815

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  AlkhalifahT. (2014). Effective wavefield extrapolation in anisotropic media: accounting for resolvable anisotropy. Geophys. Prospect.62 (5), 1089–1099. 10.1111/1365-2478.12121

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  AppelöD.PeterssonN. A. (2009). A stable finite difference method for the elastic wave equation on complex geometries with free surface. Commun. Comput. Phys.5 (1), 84–107. 10.4208/cicp.2009.v5.p84

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BaylissA.JordanK.LemesurierB. J.TurkelE. (1986). A fourth-order accurate finite-difference scheme for the computation of elastic waves. Bull. Seismol. Soc. Am.76 (4), 1115–1132. 10.1785/bssa0760041115

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  BaysalE.KosloffD. D.SherwoodJ. W. C. (1983). Reverse time migration. Geophysics48 (11), 1514–1524. 10.1190/1.1441434

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BooreD. M. (1972). “Finite-difference methods for seismic waves,” in Methods in Computational Physics. Editor BoltB. A. (New York: Academic Press), 11, 1–37.

  • Google Scholar
  • View reference in article

• 9

  BorisovD.ModrakR.GaoF.TrompJ. (2018). 3D elastic full-waveform inversion of surface waves in the presence of irregular topography using an envelope-based misfit function. Geophysics83 (1), R1–R11. 10.1190/GEO2017-0081.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  BrebbiaC. A. (1981). Progress in Boundary Element Methods, 1. Plymouth, MA: Pentech Press.

  • Google Scholar
  • View reference in article

• 11

  CerjanC.KosloffD.KosloffR.ReshefM. (1985). A nonreflecting boundary condition for discrete acoustic and elastic wave equations. Geophysics50 (4), 705–708. 10.1190/1.1441945

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  ChenZ.PanW.CaoH.SongJ.ZhangY.WangC.et al (2026). Applying VTI wave-equation envelope travel-time inversion to the walk-away VSP data in the sichuan basin of China. J. Appl. Geophys.247, 106153. 10.1016/j.jappgeo.2026.106153

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  DominguezJ. (1993). Boundary Elements in Dynamics. Southampton, UK: Computational Mechanics Publications.

  • Google Scholar
  • View reference in article

• 14

  DuveneckE.BakkerP. M. (2011). Stable P-wave modeling for reverse-time migration in tilted TI media. Geophysics76 (2), S65–S75. 10.1190/1.3533964

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  FornbergB. (1988). The pseudospectral method: accurate representation of interfaces in elastic wave calculations. Geophysics53 (5), 625–637. 10.1190/1.1442497

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  GeZ. X.FuL. Y.WuR. S. (2005). P-SV wave-field connection technique for regional wave propagation simulation. Bull. Seismol. Soc. Am.95 (4), 1375–1386. 10.1785/0120040173

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  GravesR. W. (1996). Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences. Bull. Seismol. Soc. Am.86 (4), 1091–1106. 10.1785/bssa0860041091

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  GrayS. H.MarfurtK. J. (1995). Migration from topography: improving the near-surface image. Can. J. Expl. Geophys.31 (1/2), 18–24.

  • Google Scholar
  • View reference in article

• 19

  GrechkaV.ZhangL.RectorJ. W. (2004). Shear waves in acoustic anisotropic media. Geophysics69 (2), 576–582. 10.1190/1.1707077

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  HestholmS.RuudB. (1994). 2-D finite-difference elastic wave modeling including surface topography. Geophys. Prospect.42 (5), 371–390. 10.1111/j.1365-2478.1994.tb00216.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  HestholmS.RuudB. (1998). 3-D finite-difference elastic wave modeling including surface topography. Geophysics63 (2), 613–622. 10.1190/1.1444360

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  HixonR. (1997). On Increasing the Accuracy of Maccormack Schemes for Aeroacoustic Applications. Atlanta, GA: American Institute of Aeronautics and Astronautics, 97–1586. 10.2514/6.19971586

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  HuF. Q.HussainiM. Y.MantheyJ. L. (1996). Low-dissipation and low-dispersion Runge-Kutta schemes for computational acoustics. J. Comput. Phys.124, 177–191. 10.1006/jcph.1996.0052

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  IgelH.MoraP.RiolletB. (1995). Anisotropic wave propagation through finite-difference grids. Geophysics60 (4), 1203–1216. 10.1190/1.1443849

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  JihR. S.McLaughlinK. L.DerZ. A. (1988). Free-boundary conditions of arbitrary polygonal topography in a two-dimensional explicit elastic finite-difference scheme. Geophysics53 (8), 1045–1055. 10.1190/1.1442541

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  KomatitschD.TrompJ. (1999). Introduction to the spectral element method for three-dimensional seismic wave propagation. Geophys. J. Int.139 (3), 806–822. 10.1046/j.1365-246x.1999.00967.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  KomatitschD.BarnesC.TrompJ. (2000). Simulation of anisotropic wave propagation based upon a spectral element method. Geophysics65 (4), 1251–1260. 10.1190/1.1444816

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  KristekJ.MoczoP.ArchuletaR. J. (2002). Efficient methods to simulate planar free surface in the 3D 4th-order staggered-grid finite-difference schemes. Studia Geophys. Geod.46, 355–381. 10.1023/a:1019866422821

  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  LanH. Q.ZhangZ. J. (2011). Three-dimensional wave-field simulation in heterogeneous transversely isotropic medium with irregular free surface. Bull. Seismol. Soc. Am.101 (3), 1354–1370. 10.1785/0120100194

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  LevanderA. R. (1988). Fourth-order finite-difference P-SV seismograms. Geophysics53 (11), 1425–1436. 10.1190/1.1442422

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  LiX.YaoG.NiuF.WuD. (2020). An immersed boundary method with iterative symmetric interpolation for irregular surface topography in seismic wavefield modelling. J. Geophys. Eng.17 (4), 643–660. 10.1093/jge/gxaa019

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  LiA.ZhuX.LiuT.LiaoJ.LiuX.YouJ. (2024). One-step sparse-matrix method for full wave equation depth migration. Geophysics89 (6), S405–S413. 10.1190/geo2023-0372.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  LiA.LiC.LaiF.LiaoJ.LiuT. (2025). An acoustic wave-equation depth migration method using generalized two-way phase shift operator for anisotropic media. IEEE Geoscience Remote Sens. Lett.22, 1–5. 10.1109/lgrs.2025.3549919

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  LisitsaV.VishnevskiyD. (2010). Lebedev scheme for the numerical simulation of wave propagation in 3D anisotropic elasticity. Geophys. Prospect.58 (4), 619–635. 10.1111/j.1365-2478.2009.00862.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  LiuQ. C.ZhangJ. F.GaoH. W. (2017). Reverse-time migration from rugged topography using irregular, unstructured mesh. Geophys. Prospect.65 (2), 453–466. 10.1111/1365-2478.12415

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  MacCormackR. W. (1969). The Effect of Viscosity in Hypervelocity Impact Cratering. Cincinnati, Ohio: American Institute of Aeronautics and Astronautics, 69–354.

  • Google Scholar
  • View reference in article

• 37

  MooreE. H. (1920). On the reciprocal of the general algebraic matrix, Bull. Am. Math. Soc.26, 394–395.

  • Google Scholar
  • View reference in article

• 38

  PenroseR. (1955). A generalized inverse for matrices. Math. Proc. Camb. Phil. Soc.51, 406–413. 10.1017/S0305004100030401

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  SaengerE. H.BohlenT. (2004). Finite-difference modelling of viscoelastic and anisotropic wave propagation using rotated staggered grid. Geophysics69 (2), 583–591. 10.1190/1.1707078

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  SongJ.CaoH.PanW.YangZ.LiH.LuM.et al (2023). Attenuation sensitivity kernel analysis in viscoelastic full-waveform inversion based on the generalized standard linear solid rheology. Surv. Geophys.44, 1041–1073. 10.1007/s10712-023-09767-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  SunY. C.ZhangW.ChenX. F. (2016). Seismic-wave modeling in the presence of surface topography in 2D general anisotropic media by a curvilinear grid finite-difference method. Bull. Seismol. Soc. Am.106 (3), 1036–1054. 10.1785/0120150285

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  TarantolaA. (1984). Inversion of seismic reflection data in the acoustic approximation. Geophysics49 (8), 1259–1266. 10.1190/1.1441754

  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  ThomsenL. (1986). Weak elastic anisotropy. Geophysics51 (10), 1954–1966. 10.1190/1.1442051

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  VirieuxJ. (1984). SH-wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics49 (11), 1933–1942. 10.1190/1.1441605

  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  VirieuxJ. (1986). P-SV-wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics51 (4), 889–901. 10.1190/1.1442147

  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  VirieuxJ.OpertoS. (2009). An overview of full-waveform inversion in exploration geophysics. Geophysics74 (6), WCC1–WCC26. 10.1190/1.3238367

  • CrossRef
  • Google Scholar
  • View reference in article

• 47

  WangX.LiuX. (2007). 3-D acoustic wave equation forward modeling with topography. Appl. Geophys.4 (1), 8–15. 10.1007/s11770-007-0002-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  XuS.StovasA.AlkhalifahT.MikadaH. (2020). New acoustic approximation for transversely isotropic media with a vertical symmetry axis. Geophysics85 (1), C1–C12. 10.1190/geo2019-0100.1

  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  ZhangW.ChenX. F. (2006). Traction image method for irregular free surface boundaries in finite difference seismic wave simulation. Geophys. J. Int.167 (1), 337–353. 10.1111/j.1365-246X.2006.03113.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 50

  ZhangW.ShenY. (2010). Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling. Geophysics75 (4), T141–T154. 10.1190/1.3463431

  • CrossRef
  • Google Scholar
  • View reference in article

• 51

  ZhangW.ZhangZ. G.ChenX. F. (2012). Three-dimensional elastic wave numerical modelling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids. Geophys. J. Int.190 (1), 358–378. 10.1111/j.1365-246X.2012.05472.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  ZhangZ. G.ZhangW.ChenX. F. (2014). Three-dimensional curved grid finite-difference modelling for non-planar rupture dynamics. Geophys. J. Int.199 (2), 860–879. 10.1093/gji/ggu308

  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  ZhuJ.DormanJ. (2000). Two-dimensional, three-component wave propagation in a transversely isotropic medium with arbitrary orientation: finite-element modelling. Geophysics65 (3), 934–942. 10.1190/1.1444789

  • CrossRef
  • Google Scholar
  • View reference in article

• 54

  ZhuH. J.ZhangW.ChenX. F. (2009). Two-dimensional seismic wave simulation in anisotropic media by non-staggered finite difference method. Chin. J. Geophys.52 (6), 1536–1546. 10.1190/1.1444789

  • CrossRef
  • Google Scholar
  • View reference in article