Simulate the elastic wavefields in media with an irregular surface topography based on staggered grid finite difference

Abstract

Wave equation forward modeling is a useful method to study the propagation regulation of seismic wavefields. Finite difference (FD) is one of the most extensively employed numerical approaches for computing wavefields in earthquake and exploration seismology. However, the FD approach relying on regular grids often struggles to calculate wavefields in regions featuring surface topographies. The elastic wave equation can more accurately describe the propagation of seismic wavefields in elastic media compared to the acoustic wave equation. We introduce a new FD scheme to calculate the elastic wavefields in an isotropic model with a surface topography. The novel approach can use a conventional staggered grid FD (SGFD) approach based on regular grids. A new elastic model with a horizontal surface is first obtained from the nearby surface's elastic properties and the undulating terrain elevation. We subsequently employ a topography-related strategy to eliminate the effects of surface topographies on the seismic wavefields in models with irregular surface topographies. The merits of our proposed scheme lie in its ability to stable numerically compute wavefields in models with irregular surface topographies without altering the conventional SGFD relying on regular grids. To validate the effectiveness and practicality of our method, we utilize elastic models featuring complex surface topographies. Numerical experiments demonstrate that our approach efficiently calculates elastic wavefields in isotropic media with irregular topographies based on conventional SGFD.

forward modeling, staggered grid finite difference, surface topography, elastic wavefields

1. Introduction

Wave equation forward modeling is an effective and powerful tool to study and understand the propagation of seismic wavefields in global seismology and exploration seismology. Physical analog and numerical simulation are the two main forward modeling methods. The numerical simulation is an efficient and commonly used wave equation forward modeling method at present. Among numerical techniques, three primary numerical simulation techniques have been introduced for solving wave equations: the FD approach (Kelly et al.1976), the finite element (FE) approach (Komatitsch and Vilotte 1998), and the pseudo-spectral (PS) approach (Kosloff and Baysal 1982, Fornberg 1988). The FE approach is good at simulating the wavefields in models with complicated geometries and boundaries naturally, but it comes with high computational burden and memory requirements. The PS approach uses the Fourier transform to solve the spatial partial derivative in the partial differential equation, which also contributes to significant computational expense. By contrast, the FD approach uses finite differences to approximate the spatial partial derivative, offering a more manageable computational load and memory footprint. Its straightforward implementation makes the FD approach the preferred choice for wave equation numerical simulation, particularly in wave equation inversion and imaging of exploration seismology.

Irregular topography significantly influences the propagation of seismic wavefields, necessitating specialized numerical methods for accurate calculations wavefields in models with complex terrain (Rao and Wang 2013, 2018). Common approaches include the FE approach, body-fitted grid or coordinate transformation FD approaches, and unstructured grids FD approach, which can numerically calculate seismic wavefields in models with complex surface topography. As discussed before, the classical FE approach (Komatitsch and Vilotte 1998) encounters substantial computational and memory demands, along with challenging mesh generation processes. The most straightforward and efficient FD approach for calculating the wavefields beneath irregular surface topography is the regular rectangular grid-based FD approach due to the high efficiency of the FD approach and generating a regular rectangular grid. However, conventional regular rectangular grid-based FD approach is prone to generating diffraction waves related to irregular topography, potentially compromising wavefields’ accuracy. Although increasing grid density can weaken the generation of diffracted waves, but leading to complex grid divisions and increased computational complexity. The body-fitted grid or coordinate transformation FD approach first converts the physical curvilinear grids to calculational regular rectangular grids (Hestholm and Ruud 1994, Nielsen et al.1994, 1998, Zhang and Chen 2006, Lan and Zhang 2011, Zhang et al.2012, Qu et al.2017, Liu et al.2018, Sun et al.2019), the conventional FD approach is then utilized to compute the wavefields at these calculational regular rectangular grid grids. The unstructured grids method employs advantages of the FE and FD approaches to compute wavefields in surface topography (Zhang and Liu 1999). However, it struggles with mesh extension, especially in 3D. In conclusion, the most practical and efficient approach to calculating seismic wavefields in models with irregular topographies is the body-fitted grid or coordinate transformation FD method.

There are two types of regular rectangular grid-based FD approach: collocated grid FD (CGFD) approach and SGFD approach. The CGFD method involves assigning all the wavefields variables and model parameters to the same grid, thereby rendering it a preferred technique in the body-fitted grid or coordinate transformation FD approaches. However, the CGFD approach induces an unstable checkerboard-like distribution errors (Hixon and Turkel 2000, Zhang et al.2012, Sun et al.2019). To weaken such unstable errors, a higher-order MacCormack scheme is a good choice in the body-fitted grid or coordinate transformation FD approach (Hixon and Turkel 2000, Zhang et al.2012, Sun et al.2019). The SGFD approach (Virieux 1984) sets the wavefields variables and model parameters at different grids. It is good at computing spatial derivatives in complex and highly variable heterogeneous media, effectively eliminating the unstable checkerboard-like distribution errors, and ensuring robust computation. However, the accuracy of the SGFD approach will decrease when computing seismic wavefields in models with irregular topographies, as it necessitates wavefields interpolation. To improve computational stability and accuracy of the SGFD approach, a Lebedev grid FD (Lisitsa and Vishnevskiy 2010, Qu et al.2017) scheme is commonly employed. However, the Lebedev grid FD approach requires two separate sets of staggered grids, which results in a doubling of both computational expenses and memory requirements compared to the standard SGFD approach.

There are three main issues with the coordinate transformation FD approach: first, the coordinate transformation FD approach requires the surface topography to be smooth and differentiable, so it is difficult to handle rapidly changing irregular topography. Then, the structural information will change after the coordinate transformation, which results in the decline of its forward modeling accuracy. In addition, it requires additional partial derivative calculations related to the coordinate system, which increases the calculation amount, and its calculation process is also more cumbersome. Although the accuracy of body-fitted grid FD approach is higher than that of the coordinate transformation FD approach. However, the body-fitted grid FD approach (Zhang and Chen 2006, Lan and Zhang 2011, Zhang et al.2012, Qu et al.2017, Liu et al.2018, Sun et al.2019) generates body-fitted grids by solving partial differential equations, which requires a large amount of computation. It also needs additional calculation of the partial derivative related to the coordinate system.

P- and S-waves simultaneously propagate through subsurface media, not just P-waves. Technologies related to elastic media can pinpoint oil and gas reserves, particularly in reservoirs and unique geological structures, more precisely than those associated with the acoustic approximation. Recently, elastic reverse time migration (ERTM) (Sun and McMechan 2001, Yan and Sava 2008, Du et al.2012, Rocha et al.2016, Ren et al.2017, Zhong et al.2019, Liu et al.2021, Zhong et al.2021, 2022a) and elastic full waveform inversion (EFWI) (Tarantola 1986, Brossier et al.2009, Köhn et al.2012, Ren and Liu 2016, Zhong et al.2022b) have been attracted academic attention in seismic exploration. An efficient and stable elastic wavefields numerical forward modeling approach is the cornerstone of these multi-component technologies.

The FD approach based on a regular rectangular grid offers advantages such as high efficiency, low memory needs, low computational requirements, and straightforward mesh generation. However, it can be challenging to efficiently and robustly calculate seismic wavefields in media with an irregular surface topography using the regular rectangular grid-based FD approach.

Our article presents a straightforward and stable SGFD approach to numerically calculate elastic wavefields in regions featuring surface topographies. We first fill the elastic parameters between the maximum elevation of the surface and the surface topography with the near-surface elastic parameter, thereby obtaining a new elastic model with a horizon topography. Subsequently, we use a strategy associated with the surface topography for effectively eliminating the surface topography's influence.

In the upcoming sections, we first introduce the isotropic first-order particle velocity-stress elastic wave equation, followed by the conventional SGFD method. Then, we present some schemes to mitigate the impact of surface topography on wavefields. To illustrate the effectiveness and stability of the new approach, we use some numerical examples in the description of the results. Finally, we gave some discussions and draw some conclusions.

2. Methods

Owing to its high accuracy, ability to simulate heterogeneous media, and ease of implementing absorbing boundary conditions, this research employs the isotropic particle velocity-stress elastic wave equation to forward the isotropic elastic wavefields.

2.1. Elastic wave equations in isotropic media

The 2D elastic wave equations based on particle velocity and stress (Virieux 1984), for describing wavefields propagation in an isotropic elastic model, can be expressed as follows:

$$\begin{eqnarray} \left\{ \begin{array}{@{}l@{}} \displaystyle\rho \frac{{\partial {{V}_x}}}{{\partial t}} = \frac{{\partial {{\sigma }_{xx}}}}{{\partial x}} + \frac{{\partial {{\sigma }_{xz}}}}{{\partial z}},\\ \displaystyle\rho \frac{{\partial {{V}_z}}}{{\partial t}} = \frac{{\partial {{\sigma }_{zx}}}}{{\partial x}} + \frac{{\partial {{\sigma }_{zz}}}}{{\partial z}},\\ \displaystyle\frac{{\partial {{\sigma }_{xx}}}}{{\partial t}} = (\lambda + 2\mu )\frac{{\partial {{V}_x}}}{{\partial x}} + \lambda \frac{{\partial {{V}_z}}}{{\partial z}} + {{f}_{xx}},\\ \displaystyle\frac{{\partial {{\sigma }_{zz}}}}{{\partial t}} = (\lambda + 2\mu )\frac{{\partial {{V}_z}}}{{\partial z}} + \lambda \frac{{\partial {{V}_x}}}{{\partial x}} + {{f}_{zz}},\\ \displaystyle\frac{{\partial {{\sigma }_{zx}}}}{{\partial t}} = \mu \left(\frac{{\partial {{V}_x}}}{{\partial z}} + \frac{{\partial {{V}_z}}}{{\partial x}}\right), \end{array} \right.\end{eqnarray}$$

(1)

where |${{V}_x}$| and |${{V}_z}$| are the particle velocities in horizontal and vertical directions, respectively; and |${{\sigma }_{xx}}$|⁠, |${{\sigma }_{zz}}$| and |${{\sigma }_{zx}}$| are the horizontal normal, vertical normal, and shear stress, respectively. |${{f}_{xx}}$| and |${{f}_{zz}}$| are source items. The Lamé coefficient is denoted by |$(\lambda ,\mu )$|⁠, and the density is denoted by |$\rho $|⁠. The relationship among P- and S-wave velocities, Lamé coefficient and density is expressed as such in isotropic elastic media:

$$\begin{eqnarray} {{V}_p}{\rm{ = }}\sqrt {\frac{{\lambda {\rm{ + 2}}\mu }}{\rho}} \quad {\rm{and}} \quad {{V}_s}{\rm{ = }}\sqrt {\frac{\mu }{\rho}}. \end{eqnarray}$$

2.2. Numerical scheme based SGFD to simulate elastic wavefields in an isotropic model with an irregular surface topography

The elastic wave Equations (1) are numerically solved by the SGFD approach (Virieux 1984) in this research. The discretization equations derived from the SGFD approach corresponding to the elastic wave Equations (1) are as follows:

$$\begin{eqnarray} \left\{ \begin{array}{@{}l@{}} \displaystyle{{V}_x}_{i,j}^{k + 1/2} = {{V}_x}_{i,j}^{k - 1/2} + \ \frac{{\Delta t}}{{{{\rho }_{i,j}}}}\left\{ {\frac{1}{{\Delta x}}\sum\limits_{n = 1}^N {{{C}_n}\left[ {{{\sigma }_{xx}}_{i + \left( {2n - 1} \right)/2,j}^k - {{\sigma }_{xx}}_{i - \left( {2n - 1} \right)/2,j}^k} \right] + \frac{1}{{\Delta z}}\sum\limits_{n = 1}^N {{{C}_n}\left[ {{{\sigma }_{xz}}_{i,j + \left( {2n - 1} \right)/2}^k - {{\sigma }_{xz}}_{i,j - \left( {2n - 1} \right)/2}^k} \right]} } } \right\},\\ \displaystyle{{V}_z}_{i + 1/2,j + 1/2}^{k + 1/2} = {{V}_z}_{i + 1/2,j + 1/2}^{k - 1/2} + \frac{{\Delta t}}{{{{\rho }_{i + 1/2,j + 1/2}}}}\left\{ {\frac{1}{{\Delta x}}\sum\limits_{n = 1}^N {{{C}_n}\left[ {{{\sigma }_{xz}}_{i + n,j + 1/2}^k - {{\sigma }_{xz}}_{i - \left( {n - 1} \right),j + 1/2}^k} \right] + \frac{1}{{\Delta z}}\sum\limits_{n = 1}^N {{{C}_n}\left[ {{{\sigma }_{zz}}_{i + 1/2,j + n}^k - {{\sigma }_{zz}}_{i + 1/2,j - \left( {n - 1} \right)}^k} \right]} } } \right\},\\ \displaystyle{{\sigma }_{xx}}_{i + 1/2,j}^{k + 1} = {{\sigma }_{xx}}_{i + 1/2,j}^k + \ \Delta t\left\{ {\frac{{{{{\left( {\lambda + 2\mu } \right)}}_{i + 1/2,j}}}}{{\Delta x}}\sum\limits_{n = 1}^N {{{C}_n}\left[ {{{V}_x}_{i + n,j}^{k + 1/2} - {{V}_x}_{i - \left( {n - 1} \right),j}^{k + 1/2}} \right] + \frac{{{{\lambda }_{i + 1/2,j}}}}{{\Delta z}}\sum\limits_{n = 1}^N {{{C}_n}\left[ {{{V}_z}_{i + 1/2,j + (2n - 1)/2}^{k + 1/2} - {{V}_z}_{i + 1/2,j - \left( {2n - 1} \right)/2}^{k + 1/2}} \right]} } } \right\} + {{f}_{xx}},\\ \displaystyle{{\sigma }_{zz}}_{i + 1/2,j}^{k + 1} = {{\sigma }_{zz}}_{i + 1/2,j}^k + \Delta t\left\{ {\frac{{{{{\left( {\lambda + 2\mu } \right)}}_{i + 1/2,j}}}}{{\Delta z}}\sum\limits_{n = 1}^N {{{C}_n}\left[ {{{V}_z}_{i + 1/2,j + \left( {2n - 1} \right)/2}^{k + 1/2} - {{V}_z}_{i + 1/2,j - \left( {2n - 1} \right)/2}^{k + 1/2}} \right] + \frac{{{{\lambda }_{i + 1/2,j}}}}{{\Delta x}}\sum\limits_{n = 1}^N {{{C}_n}\left[ {{{V}_x}_{i + n,j}^{k + 1/2} - {{V}_x}_{i - \left( {n - 1} \right),j}^{k + 1/2}} \right]} } } \right\} + {{f}_{zz}},\\ \displaystyle{{\sigma }_{xz}}_{i,j + 1/2}^{k + 1} = {{\sigma }_{xz}}_{i,j + 1/2}^k + \Delta t\left\{ {\frac{{{{\mu }_{i,j + 1/2}}}}{{\Delta z}}\sum\limits_{n = 1}^N {{{C}_n}\left[ {{{V}_x}_{i,j + n}^{k + 1/2} - {{V}_x}_{i,j - \left( {n - 1} \right)}^{k + 1/2}} \right]} + \frac{{{{\mu }_{i,j + 1/2}}}}{{\Delta x}}\sum\limits_{n = 1}^N {{{C}_n}\left[ {{{V}_z}_{i + \left( {2n - 1} \right)/2,j + 1/2}^{k + 1/2} - {{V}_z}_{i - \left( {2n - 1} \right)/2,j + 1/2}^{k + 1/2}} \right]} } \right\}, \end{array} \right. \end{eqnarray}$$

(2)

where k symbolizes a discrete time, i represents the grid position for x-axis discretization, and j denotes the grid position for z-axis discretization; |$\Delta t$| stands for the time step, |$\Delta x$| indicates the grid interval for x-axis, and |$\Delta z$| signifies the grid interval for z-axis. The FD coefficient |${{C}_n}$|⁠, often calculated by Taylor expansion, can also be obtained using an optimization approach (Liu 2014).

To mitigate boundary reflections, the perfectly matched layer (PML), as introduced by Berenger (1994), is employed in this research.

A filter associated with the surface topography is presented to continue the use conventional SGFD approach to calculate the wavefields in models with complex topographies. The filter associated with the surface topography has the following expression:

$$\begin{eqnarray} f(x,z) = \left\{ {\begin{array}{@{}*{2}{c}@{}} 0,&\quad if(z < \textit{topography}\, (x,z)),\\ 1,&\quad if(z \ge \textit{topography}\, (x,z)), \end{array}} \right.\end{eqnarray}$$

(3)

where the proposed topography-related filter is |$f(x,z)$|⁠, and the surface topography function is |$topography\, (x,z)$|⁠.

To obtain stable wavefields by the proposed method, we first construct a new model featuring a horizontal surface. The elastic parameters between the surface's highest elevation and topography are filled with the near-surface elastic parameters from the original model containing a complex surface topography. The new model, featuring a horizontal surface, is then obtained. Following that, the wavefields can be numerically computed by the conventional regular rectangular grids-based SGFD approach in the newly constructed horizontal surface model.

The final wavefields in original model, featuring complex surface topography, can be simulated by following equation:

$$\begin{eqnarray} U(x,z,t) = u(x,z,t)f(x,z), \end{eqnarray}$$

(4)

where |$U(x,z,t)$| is the actual wavefields that needs to be calculated in the original surface topography model featuring a surface topography, and |$u(x,z,t)$| is the calculated wavefields in the newly constructed model featuring a horizontal surface.

In this scheme, the actual wavefields |$U(x,z,t)$| can be easily computed utilizing the conventional SGFD approach. The proposed scheme almost requires no modifications to the conventional SGFD code, which is often used in the horizontal surface model. Moreover, the computational burden and memory requirements have essentially remained unchanged.

The flowchart of the proposed method is illustrated in Fig. 1. The comprehensive workflow contains the following steps:

Acquire model information including the elastic and terrain parameters. And generate the new model featuring a horizontal surface utilizing the near-surface elastic parameters and surface topography information.
Calculate the elastic wavefields using the discrete Equation (2) in the newly filled model via the SGFD method, and apply the PML boundary condition to eliminate boundary reflections.
Establish a topography-related filter (Equation 3) based on the surface topography parameters.
Calculate the elastic wavefields in models featuring complex surface topography using Equation (4).

Figure 1.

The flowchart of elastic wave equation forward modeling in models with complex surface topographies based on regular rectangular grids.

Open in new tab Download slide

Figure 2 displays a layered model with a topography to illustrate the proposed method to obtain a new horizontal surface model for subsequent wavefields calculations. The procedure is straightforward; fill in the velocity and density parameters near the surface between the maximum and minimum elevations of the undulating surface to build a new horizontal surface model. In this new model, conventional SGFD methods based on regular rectangular grids can be employed for numerical calculate the wavefields.

Figure 2.

The new horizontal surface obtained from the original model's surface topography. The surface topography of original model is denoted by the red line.

Open in new tab Download slide

3. Numerical examples

This section employs some elastic models to assess the effectiveness and rationality of our proposed method in this research. The conventional SGFD method, with second-order time accuracy and twelfth-order space accuracy, is utilized to numerically solve the Equations (1).

3.1. The isotropic homogeneous models with different surface topographies

A homogeneous elastic model is first used to simply and clearly demonstrate the validity of our proposed schemes. The P- and S-wave velocities of the homogeneous elastic model are 3000 and 1000 m/s, and the density is 1000 kg/m³. The maximum number of grids in both horizontal and vertical directions is 201. The grid interval in both horizontal and vertical directions is 5 m. We use a Ricker wavelet with a 30 Hz dominant frequency and a temporal interval of 0.5 ms as the source.

Two distinct surface topographies are employed to validate the effectiveness and rationality of our proposed approach: one features a smooth surface topography, where wavefields may be accurately computed utilizing the coordinate transformation FD method; the other has an uneven surface with some irregular points, for which the coordinate transformation FD method struggles to calculate high-precision wavefields in such a model.

We first acquire the new horizontal surface elastic model using the presented approach. Subsequently, a conventional SGFD method is employed to compute the wavefields in the new elastic model. For comparison, we also calculate the elastic wavefields in the horizontal surface elastic model.

Figures 3 and 4 display the elastic wavefields in homogeneous isotropic models with varying surface topographies at 200 and 250 ms, respectively. All forward wavefields in these figures have no boundary reflections, as we employed a PML to mitigate reflections during the modeling process. The forward wavefields have no artificial reflections associated with the surface topographies. Noticeably, the elastic wavefields in the homogeneous isotropic models with smooth (Figs. 3b and e, and 4b and e) and unsmooth (Figs. 3c and f, and 4c and f) surface topographies align well with those in the horizontal surface model (Figs. 3a and d, and 4a and d), except in areas of irregular surface topography. This suggests that the proposed method efficiently calculates the elastic wavefields in regions with irregular surface topography.

Figure 3.

Elastic wavefields at 200 ms in different homogeneous isotropic models. The top rows (a, b, and c) display horizontal component particle velocity wavefields, and the bottom rows (d, e, and f) exhibit vertical component particle velocity wavefields. Columns (a) and (d) represent elastic wavefields in a horizontal surface homogeneous model. Columns (b) and (e) correspond to elastic wavefields in a smooth surface homogeneous model, calculated using the suggested approach. Columns (c) and (f) show elastic wavefields in an unsmooth surface homogeneous model, also computed by the proposed method. The red line indicates the surface topography line, and the green arrow points to the unsmooth singularities.

Open in new tab Download slide

Figure 4.

Elastic wavefields at 250 ms in homogeneous isotropic models. The top rows (a, b, and c) display horizontal component particle velocity wavefields, and the bottom rows (d, e, and f) exhibit vertical component particle velocity wavefields. Columns (a) and (d) represent elastic wavefields in a horizontal surface homogeneous model. Columns (b) and (e) correspond to elastic wavefields in a smooth surface homogeneous model, calculated using the suggested approach. Columns (c) and (f) show elastic wavefields in an unsmooth surface homogeneous model, also computed by the proposed method.

Open in new tab Download slide

3.2. The isotropic graben models with different surface topographies

We further employ some graben models to show the efficacy of the proposed approach. We incorporate two different categories of surface topographies into the graben model. The first category features a smooth and differentiable surface topography (Fig. 5). We use the coordinate transformation FD technique. The second category of surface topography features irregularities and singularities (Fig. 6), causing the coordinate transformation FD method to struggle with computing high-precision wavefields in the corresponding model.

Figure 5.

A graben model featuring a smooth surface topography: (a) P-wave velocity model, (b) S-wave velocity model, (c) density model, and (d) the surface topography of the model.

Open in new tab Download slide

Figure 6.

A graben model with an unsmooth surface topography: (a) P-wave velocity model, (b) S-wave velocity model, (c) density model, and (d) the surface topography of the model.

Open in new tab Download slide

Our proposed method yields the new horizontal surface elastic model, then the wavefields are calculated in this new model. For comparison, we also generate reference wavefields by calculating the wavefields in the horizontal surface elastic model. We calculate the acoustic wavefields to show the different characteristics of acoustic and elastic wavefields.

Figure 7 displays the different forward acoustic wavefields in different graben models; whereas the corresponding elastic wavefields are shown in Fig. 8. We can see the elastic wavefields are more complex than acoustic wavefields from Figs. 7 and 8. Because P- and S-waves are mixed in elastic wavefields, while there are only P-waves in acoustic wavefields. We can observe that the propagation of acoustic and elastic wavefields changes with the surface topographies. The acoustic and elastic wavefields do not exhibit any false diffraction waves associated with topographies.

Figure 7.

Acoustic wavefields at 600 ms in different graben models by the proposed scheme. The top rows (a, b, and c) display horizontal component particle velocity wavefields, and the bottom rows (d, e, and f) exhibit vertical component particle velocity wavefields. Columns (a) and (d) represent acoustic wavefields in the horizontal surface model. Columns (b) and (e) correspond to acoustic wavefields in the smooth surface model, calculated using the suggested approach. Columns (c) and (f) show acoustic wavefields in the unsmooth surface model, also computed by the proposed method.

Open in new tab Download slide

Figure 8.

Elastic wavefields at 600 ms in different graben models by the proposed scheme. The top rows (a, b, and c) display horizontal component particle velocity wavefields, and the bottom rows (d, e, and f) exhibit vertical component particle velocity wavefields. Columns (a) and (d) represent elastic wavefields in the horizontal surface model. Columns (b) and (e) correspond to elastic wavefields in the smooth surface model, calculated using the suggested approach. Columns (c) and (f) show elastic wavefields in the unsmooth surface model, also computed by the proposed method.

Open in new tab Download slide

Figures 9 and 10 present the corresponding acoustic and elastic multi-component synthetic seismic records. In Figs. 9 and 10, we can observe that the propagation of multi-component synthetic seismic records changes with the surface topographies. In addition, there are no any false diffraction waves associated with topographies in multi-component synthetic seismic records.

Figure 9.

Different acoustic synthetic seismic records in the graben models with different surface topographies. The horizontal component synthetic seismic records are shown in top rows (a, b, and c), and the vertical component synthetic seismic records are shown in bottom rows (d, e, and f). Acoustic synthetic seismic records in the horizontal surface model are shown in (a) and (d). Acoustic synthetic seismic records in the smooth surface model, which are calculated by the proposed scheme, are shown in (b) and (e); and acoustic synthetic seismic records in the unsmooth surface model, which are calculated by the proposed scheme, are shown in (c) and (f).

Open in new tab Download slide

Figure 10.

Different elastic multi-component synthetic seismic records in the graben models with different surface topographies. The horizontal component synthetic seismic records are shown in top rows (a, b, and c), and the vertical component synthetic seismic records are shown in bottom rows (d, e, and f). Multi-component synthetic seismic records in the horizontal surface model are shown in (a) and (d). Multi-component synthetic seismic records in the smooth surface model, which are calculated by the proposed scheme, are shown in (b) and (e); and multi-component synthetic seismic records in the unsmooth surface model, which are calculated by the proposed scheme, are shown in (c) and (f).

Open in new tab Download slide

Owing to alterations in wavefields propagation caused by surface topography, reflection events exhibit a non-hyperbolic pattern (Fig. 9b, c, e and f), differing from models with horizontal topography (Fig. 9a and d). These events in models with varying surface topographies change accordingly. Elastic multi-component recording also presents a similar situation (Fig. 10). We computed the misfit between elastic waveforms in the horizontal topography model and those in the complex topography models. Figure 11 displays the extracted traces at 1000 m from Fig. 10. Since the source and receiver in the complex topography models are at the same 1000 m location, the waveforms in the horizontal topography model align well with those in the other models in Fig. 10. The discrepancies in waveforms between the horizontal topography model and the complex topography models are equal to zero.

Figure 11.

The extracted elastic multi-component synthetic seismic records at 1000 m in Fig. 10. The top rows (a and b) are the horizontal component seismic records, and the bottom rows (c and d) are the vertical component seismic records. The red solid lines display waveforms in the horizontal surface model, black circled lines represent waveforms in the smooth surface model, blue circled lines show waveforms in the unsmooth surface model, and green lines indicate the difference waveforms between the waveforms in the horizontal surface model and those in models with different complex surface topographies. ‘HS’, ‘SST’, ‘DHS’, ‘USST’, and ‘DHU’ represent waveforms in the horizontal surface model, the smooth surface model, the difference between waveforms in the horizontal and smooth surface topography models, the waveforms in the unsmooth surface model, and the difference between waveforms in the horizontal and unsmooth surface topography models, respectively.

Open in new tab Download slide

This test shows the robustness and effectiveness of the proposed scheme in numerically calculating elastic wavefields within isotropic media with complex surface topographies.

3.3. The modified isotropic Hess model with different surface topographies

Finally, a more complex modified Hess model is employed to assess the effectiveness and feasibility of the proposed approach.

Figure 12 displays the modified Hess model with a smooth topography, and Fig. 13 illustrates the model with an unsmooth topography. The proposed method yields a new horizontal topography variant of the modified Hess model. The source wavelet's peak frequency is 15 Hz, and 401 multi-component geophones are positioned at the model's surface. The source position is |$( {x = 2000\ {\rm m},\ z = 0\ {\rm m}} )$|⁠.

Figure 12.

A modified Hess model featuring a smooth surface topography: (a) P-wave velocity model, (b) S-wave velocity model, (c) density model, and (d) the surface topography of the model.

Open in new tab Download slide

Figure 13.

A modified Hess model featuring an unsmooth surface topography: (a) P-wave velocity model, (b) S-wave velocity model, (c) density model, and (d) the surface topography of the model.

Open in new tab Download slide

Figure 14 displays the different forward elastic wavefields at 1100 ms. Figure 15 shows the elastic multi-component synthetic seismic records in different models. There are no false diffraction waves in Fig. 14. Note the change in travel time of reflection events with the surface topographies (Fig. 15b, c, e and f), in comparison to Fig. 15a and d. However, no surface topography-related artificial reflections are presented in the seismic records in the models with different topographies. We extracted the traces at 2000 m in Fig. 15 to show more details about the wavefields. Theoretically, records in the different complex topography models should be the same at 2000 m of the surface, as the source and receiver locations are identical in these models. We computed the misfit between the waveforms in the horizontal topography model and those in the model with different complex topographies. In Fig. 16, we can observe that waveforms in the model with horizontal topography are apparently in good agreement with waveforms in the model with different complex topographies. The differences in waveforms between those in the horizontal surface model and those in the models with different complex topographies are also equal to zero.

Figure 14.

Elastic wavefields at 1100 ms in the modified Hess models with different surface topographies. The top rows (a, b, and c) display horizontal component particle velocity wavefields, while the bottom rows (d, e, and f) exhibit vertical component particle velocity wavefields. Columns (a) and (d) represent elastic wavefields in the horizontal surface model. Columns (b) and (e) correspond to elastic wavefields in the smooth surface model, calculated using the suggested approach. Columns (c) and (f) show elastic wavefields in the unsmooth surface model, also computed by the proposed method.

Open in new tab Download slide

Figure 15.

Elastic multi-component synthetic seismic records in the modified Hess models with different surface topographies. The top rows (a, b, and c) are the horizontal component seismic records, and the bottom rows (d, e, and f) are the vertical component seismic records. (a and d) The seismic records in the horizontal surface model. (b and e) The seismic records in the smooth surface model, calculated using the suggested approach. (c and f) The seismic records in the unsmooth surface model, also computed by the proposed method.

Open in new tab Download slide

Figure 16.

The extracted elastic multi-component synthetic seismic records at 1000 m in Fig. 15. The top rows (a and b) are the horizontal component seismic records, and the bottom rows (c and d) are the vertical component seismic records. The red solid lines display waveforms in the horizontal surface model, black circled lines represent waveforms in the smooth surface model, blue circled lines show waveforms in the unsmooth surface model, and green lines indicate the difference waveforms between the waveforms in the horizontal surface model and those in models with different complex surface topographies. ‘HS’, ‘SST’, ‘DHS’, ‘USST’, and ‘DHU’ represent waveforms in the horizontal surface model, the smooth surface model, the difference between waveforms in the horizontal and smooth surface topography models, the waveforms in the unsmooth surface model, and the difference between waveforms in the horizontal and unsmooth surface topography models, respectively.

Open in new tab Download slide

The proposed numerical scheme based on SGFD to simulate wavefields in an isotropic model with an irregular surface topography can further be used to wave separation ERTM (Du et al.2012, Zhong et al.2019, Zhong et al.2021). Then, combing with the numerical scheme based on SGFD in media with an irregular surface topography, we use the wave separation ERTM based on decoupled equations to migrate the isotropic multi-component seismic data in the models with different surface topographies for further demonstration of the advantages of the newly proposed method. The final migration results are shown in Fig. 17. We cannot see false artificial reflections related to the irregular surface topography in the migration results of isotropic ERTM from Fig. 17. This implies that the topography isotropic ERTM based our proposed numerical scheme based on SGFD is effective in eliminating the impact of irregular surface topography on the images of ERTM.

Figure 17.

Migration images of the elastic multi-component seismic data in modified Hess elastic models with different irregular surface topographies obtained by the wave separation ERTM. The PP images in modified Hess elastic model with a smooth irregular topography are shown in (a); the PS images in modified Hess elastic model with a smooth irregular topography are shown in (b); the PP images in modified Hess elastic model with an unsmooth irregular topography are shown in (c); and the PS images in modified Hess elastic model with an unsmooth irregular topography are shown in (d).

Open in new tab Download slide

This complex modified Hess model effectively illustrates the credibility and rationality of our proposed approach for calculating wavefields in complex elastic structures with different complex topographies.

4. Conclusions

We have introduced a novel and efficient forward approach to compute the elastic wavefields in isotropic media with surface topography. This proposed method requires only a minimal addition of grids in the vertical direction for complex topography models. The conventional SGFD method remains unchanged, enabling efficient computation in areas with irregular surface topography. Numerical experiments illustrate that the presented approach can efficiently calculate elastic wavefields in isotropic media with surface topographies.

Conflict of interest statement. None declared.

Funding

This research was funded by National Natural Science Foundation of China (grant nos. 42204125 and 12304357), Hubei Provincial Natural Science Foundation (grant nos. 2024AFB1043, 2023AFB891, and 2023AFB352), Doctoral Research Startup Fund Project of Hubei University of Automotive Technology (grant nos. BK202308 and BK202210), Open Fund of Key Laboratory of Exploration Technologies for Oil and Gas Resources (Yangtze University), Ministry of Education (grant no. K2023-10), Chunhui Project Foundation of the Education Department of China (grant no. HZKY20220337), Excellent Youth Fund of Hubei University of Automotive Technology (grant no. 2023YQ01), Shandong Provincial Natural Science Foundation (grant no. ZR2023MD116), and Qingdao Natural Science Foundation (grant no. 23-2-1-210-zyyd-jch). Shiyan Key Laboratory of Quantum Information and Precision Optics (Hubei University of Automotive Technology) (grant no. SYZDK12024A02), the program of outstanding young and middle-aged science and technology innovation team in Hubei Province (grant no. T2023012).

Data availability

Data associated with this research are available and can be obtained by contacting the corresponding author upon reasonable request.

References

Berenger

A perfectly matched layer for the absorption of electromagnetic wave

J Comput Phys

1994

;

114

185

–

200

Month:	Total Views:
July 2024	68
August 2024	43
September 2024	116
October 2024	72
November 2024	64
December 2024	30
January 2025	24
February 2025	20
March 2025	119
April 2025	50

Article Contents

Simulate the elastic wavefields in media with an irregular surface topography based on staggered grid finite difference

Abstract

1. Introduction

2. Methods

2.1. Elastic wave equations in isotropic media

2.2. Numerical scheme based SGFD to simulate elastic wavefields in an isotropic model with an irregular surface topography

3. Numerical examples

3.1. The isotropic homogeneous models with different surface topographies

3.2. The isotropic graben models with different surface topographies

3.3. The modified isotropic Hess model with different surface topographies

4. Conclusions

Funding

Data availability

References

Citations

Views

Altmetric

Email alerts

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only