TY - JOUR
AU1 - Li,, Xiang
AU2 - Yao,, Gang
AU3 - Niu,, Fenglin
AU4 - Wu,, Di
AB - Abstract The irregular free surface topography has a significant impact on simulations of seismic wave propagation. Therefore, an accurate representation of the irregular free surface is required for an accurate wavefield simulation. We propose an immersed boundary method used in fluid dynamics calculation to simulate acoustic waves with finite-difference in media with irregular surfaces. First, we set the number of ghost layers to half the length of the finite-difference stencil. Then, we define mirror points by orthogonally projecting the ghost points to fractional points below the free surface. We calculate the wavefield at these mirror points using an iterative symmetric interpolation method. Finally, we set the wavefield at the ghost points to the negative value of the wavefield of their corresponding mirror points. The proposed iterative symmetric interpolation method allows computing the wavefield at the mirror points more accurately and stably than the conventional immersed boundary methods. Numerical examples validate the accuracy and stability of this method in seismic forward modelling with strongly varying topography. immersed boundary method, iterative symmetric interpolation, irregular surface topography, rectangular meshes, finite-difference modelling 1. Introduction The irregular free surface has a significant impact on seismic data processing, including seismic migration (Beasley & Lynn 1989; Reshef 1991; He et al. 2002) and full-waveform inversion (FWI) (Tarantola 1984; Pratt 1999). If the Earth model is discretised with rectangular grids, the irregular surface will be truncated to the nearest integer points to form a staircase-like surface, which produces artificial diffractions at the edges of staircases during forward modelling as well as travel-time errors in surface-related events (Lombard & Piraux 2004). Consequently, the free surface with strong elevation variation can affect the migration result even after static correction, which shifts the recorded data to a reference datum (Shtivelman & Canning 1988). The artificial diffractions and the travel-time errors can also produce large errors in FWI (Bleibinhaus & Rondenay 2009). Studies showed that ignorance of topography larger than half the minimum wavelength could lead to significant artifacts in the inversion results (Nuber et al. 2016; Borisov et al. 2018). Consequently, successful FWI applications, especially for land seismic data, need efficient and accurate wavefield modelling algorithms dealing with irregular surface topography. Presently the main wavefield simulation methods include finite-difference methods (FDMs) (Alterman & Karal 1968; Virieux 1986; Levander 1988; Yao et al. 2016; Mulder & Huiskes 2017; Cao & Chen 2018; 2018b), pseudo-spectral methods (PSM) (Kosloff & Baysal 1982), finite-element methods (FEMs) (Marfurt 1984; Liu et al. 2014; Meng & Fu 2017) and spectral-element methods (SEMs) (Kelly et al. 1976; Komatitsch & Vilotte 1998; Oral et al. 2019). FEMs use the triangular/tetrahedral mesh to fit the irregular surface easily while SEMs employ the quadrilateral/hexahedral mesh to match the irregular surface effectively. In addition, the FEM/SEM methods can easily achieve a spatially variable grid from fine to coarse, which is particularly useful for simulating the low-velocity layer. Although FEM/SEM can use flexible meshes to deal with irregular surface topography, they are far more expensive than FDM/PSM. However, FDM/PSM with conventional rectangular grids generates artificially high-frequency scattering energy in cases where irregular surface topography exists. The simplest way to mitigate the scattering energy is to use a fine mesh for a better approximation to the surface (Hayashi et al. 2001). Numerical experiments have shown that it needs at least 60 sampling points per minimum wavelength to accurately simulate the seismic wavefield with an irregular free surface (Bohlen & Saenger 2006; Bleibinhaus & Rondenay 2009). However, a fine mesh increases the computational cost significantly. There are several other solutions to deal with the irregular surface. Vertically variable grids can reduce the artificial scattering energy as well as maintain the efficiency (Jastram & Tessmer 2010; Li et al. 2016; Wang et al. 2019). This method basically uses a fine grid in the low-velocity layer near the surface but a coarse grid elsewhere. Wavefield interpolation is used to connect the coarse mesh to the fine one in transition zones (Hayashi et al. 2001; Jastram & Tessmer 2010). The mapping methods transform a smoothly irregular surface, which is differentiable, to a flat one through a coordinate transform, which could link the physical domain to the computational domain. The rest part of the mesh is recalibrated along the depth direction based on the transformed surface. Afterward, the normal FDM/PSM methods carry out the numerical wavefield simulations in the computational domain (Tessmer et al. 1992; Nielsen et al. 1994). Unlike the mapping methods, the boundary-conforming grids (body-fitted grids) method could adapt more complex surfaces and interfaces that are not smooth or differentiable. In this method, the curve coordinate conforms to the shape of the surfaces and the interfaces in the physical domain. This method has been applied to simulate the elastic wavefield (Zhang & Chen 2006; Zhang et al. 2012; Sun et al. 2016, 2017) and the acoustic wavefield (Rao & Wang 2013, 2018) with irregular interfaces. The Earth's surface is not only irregular but is also a free surface that reflects seismic waves. There are various methods to implement the free surface condition in seismic wavefield modelling. The simplest method is the conventional vacuum modelling method (Kelly, 1976; Virieux 1986; Muir et al. 1992), which sets the elastic parameters to a small number approximating to air parameters above the free surface, as well as a small density value to avoid instability due to division by zero in the air layer. However, when the dip angle between the boundary and the mesh becomes large, the accuracy of the vacuum method will decrease rapidly with a coarse grid. The numerical test showed that it requires 60 grid points per minimum wavelength to achieve high accuracy for all dip angles in seismic wavefield modelling (Zahradník et al. 1993; Graves 1996; Bohlen & Saenger 2006). It also increases the computational cost significantly, and therefore makes imaging processes, such as RTM and FWI, impractical. Oprsal & Zahradník (1999) extended the vacuum modelling method with non-uniform grids, in which a fine grid is used around the surface and a coarse grid is employed elsewhere. Antisymmetrically mirroring the stress components is another simple and effective method to implement the free surface boundary condition (Levander 1988). Basically, this method sets the stress components at the points above the free surface to the negative value of the stress at their mirror points and also sets the stress at the free surface to zero. The number of the involved points above the free surface in this process is half the finite-difference stencil length. Consequently, the free surface imaging condition is implemented. Zhang & Chen (2006) further developed this concept into the traction image method for elastic wavefield modelling with an irregular free surface. First, they transformed a physical grid with an irregular free surface into a rectangular grid in the computational domain, and then applied the traction image method for implementing the free surface condition. The immersed boundary method is an effective means to deal with complex interface problems. This method was first proposed for mechanical analysis in fluid dynamics for simulating cardiac dynamics and associated blood flow by Peskin (1972). So far, the immersed boundary method has been mainly employed in conjunction with Cartesian-grid Navier–Stokes solvers for flow modelling. The governing differential equations are discretised in a Cartesian coordinate system while the interface boundary on fractional grids is immersed in rectangular meshes. Different boundary conditions can be implemented by changing the field value on the neighboring integer grid. Tseng & Ferziger (2003) proposed an effective ghost point immersed boundary method for simulating turbulence in a complex geometry, where the boundary condition is enforced through the ghost cells outside the boundary. Because the grid calculation associated with the irregular surface needs to be found only once for the whole calculation, it is capable of maintaiingn the efficiency and accuracy of FDMs. The immersed boundary method was then extended for applications in other types of fields, including electric, magnetic and seismic wavefields (Lombard & Piraux 2004; Lombard et al.2008; Zhao 2010; Almuhaidib & Toksoz 2015; Gao et al. 2015; Hu 2016). Zhao (2010) used the ray casting immersed boundary method for modelling curved perfectly electric conducting walls based on Maxwell's equations. Seismic wavefield modelling with the immersed boundary method has also been performed on conventional rectangular meshes. Lombard & Piraux (2004); Lombard et al. (2008) used fictitious values built from neighboring grid nodes in a vacuum and injected them into the finite-difference calculation to modelling wave propagation with an irregular surface. Gao et al. (2015) extended the immersed boundary with staggered-grid spatial discretisation for seismic wave simulation by extrapolating the values across the free surface. One important step of the immersed boundary method is how to compute the wavefield at the ghost points according to the wavefield below the free surface. Almuhaidib & Toksoz (2015) used both interpolation and extrapolation to obtain the wavefield at the mirror points. If a mirror point is surrounded by the grid points below the surface (figure 2a of Almuhaidib & Toksoz 2015), interpolation is used to calculate the wavefield at the mirror point. Otherwise, extrapolation is used (figure 2b of Almuhaidib & Toksoz 2015). Hu (2016) used interpolation to compute the wavefield at the fictitious points, which are on a straight line with the ghost point orthogonal to the surface, and then iteratively extrapolated the wavefield from the fictitious points to the ghost points. However, both methods involve extrapolation calculations in the modelling process; therefore, it is easy for it to become unstable as well as compromise accuracy in the finite-difference modelling process. In this paper, we further improve the immersed boundary method with iterative symmetric interpolation to achieve seismic wavefield modelling with an irregular free surface topography. First, the number of ghost layers is set to half the length of the finite-difference stencil. Second, the mirror points of ghost points are found by orthogonally projecting the ghost points to fractional points below the free surface. Third, the wavefield at these mirror points is calculated by using an iterative symmetric interpolation method. Finally, the wavefield at these ghost points is set to the negative value of the wavefield at the corresponding mirror points. As extrapolation is not used in the computation of the wavefield at ghost points, our immersed boundary method is stable and accurate. Our aim is for the simulation of land seismic data on the free surface; thus, we convert the pressure wavefield into the particle vibration velocity. Numerical examples validate the accuracy and stability of this method. 2. Methodology The acoustic wave equation can be expressed as $$\begin{equation} \frac{1}{{v{{({{\bf x}})}^2}}}\frac{{{\partial ^2}p({{\bf x}},t)}}{{\partial {t^2}}} - \rho ({{\bf x}})\nabla \cdot \frac{1}{{\rho ({{\bf x}})}}\nabla p({{\bf x}},t) = \delta ({{\bf x}} - {{{\bf x}}_s})s(t), \end{equation}$$(1) where |$p({{\bf x}},t)$| represents the pressure field at the location |${{\bf x}}$| and the time t, |$s(t)$| denotes the source signature, |$v({{\bf x}})$| and |$\rho ({{\bf x}})$| are the velocity and density respectively, |$\delta ({{\bf x}} - {{{\bf x}}_s})$| represents the Dirichlet delta function for injecting the point source at the source location |${{{\bf x}}_s}$|⁠. Note that both velocity and density can vary spatially in this study. In the case of seismic wavefields propagating in the Earth, which is a half-space, the solution of equation 1 is uniquely determined by the initial condition, $$\begin{equation} p({{\bf x}},t) = 0,\,\,\,\,t < 0, \end{equation}$$(2) and the free surface boundary condition, $$\begin{equation} p({{\bf x}},t) = 0,\,\,\,{{\bf x}} \in \sum , \end{equation}$$(3) where |$\sum $| represents the free surface. The free surface separating the air and the rest of the Earth has a reflection coefficient very close to −1. In this paper, we use the acoustic wave equation 1 to simulate seismic wavefields because acoustic wave equations still prevail for seismic migration and FWI in the industry by considering two factors: low computational cost and acceptable accuracy. Acoustic wavefield modelling has one to two orders of magnitude lower computational cost than elastic wavefield modelling mainly due to the slow S-wave velocity, especially for uncompacted sediments. Acoustic wave equations can produce accurate travel-time of seismic P-wave data, which guarantees the validity of the migration algorithms based on acoustic wave equations. Also, acoustic wave equations generate relatively accurate waveform of the transmitted seismic P-wave data, which drives FWI for background updates (Warner et al. 2013; Yao et al. 2019). 2.1. Modelling the free surface with the immersed boundary method The immersed boundary method is an effective way to process an irregular free surface, and the method was first used for the mechanical analysis on the irregular boundary in the fluid dynamics field (Peskin 1972). The ‘immersion’ mainly refers to immersing the irregular boundary into a conventional Cartesian grid. The modelling process is performed on a Cartesian rectangular mesh. Therefore, with an immersed boundary method, the irregular boundary is located on fractional grids while the computational variables (e.g. the pressure wavefield) are on integer grids. The immersed boundary method for simulating an irregular free surface in seismic wavefield modelling includes several steps. First, the number of ghost layers above the irregular boundary is set to half the length of the finite-difference stencil. For instance, the number of ghost layers should be set to four for an eighth-order finite-difference modelling scheme. The computation points above the free surface are called ghost points, which are illustrated in figure 1. Second, find the mirror points of all ghost points through projection orthogonally to the surface. The middle point between a pair of mirror and ghost points is located exactly on the free surface. The middle point is referred to as an intercept point. One easy way to find an intercept point is to numerically search the point that has the shortest distance point from the ghost point to the surface boundary (Mittal et al. 2008). Since the location of the intercept points needs to be found only once during the whole calculation process, we can sample the irregular boundary densely and use the closest sampling point to the ghost point as the intercept point. The coordinates of the mirror points then can be easily obtained. Finally, the wavefield value at the ghost point is set as the negative value of the corresponding mirror point before each time step of wavefield modelling. In the case of a flat free surface, the immersed boundary method described is equivalent to the classic mirror method. Figure 1. Open in new tabDownload slide Schematic diagram of the immersed boundary method in the 2D space domain. The black solid line represents the irregular free surface boundary. The unfilled circle represents the ghost point on an integer grid above the boundary while the filled circle indicates the mirror point on a fractional grid below the boundary. The filled squares represent the points on the integer grids below the surface while the unfilled squares indicate the points on the integer grids above the surface. The straight line across a pair of the mirror and ghost points is orthogonal to the free surface. In the immersed boundary method, the wavefield value of a ghost point is the negative value of the wavefield at the corresponding mirror point. Figure 1. Open in new tabDownload slide Schematic diagram of the immersed boundary method in the 2D space domain. The black solid line represents the irregular free surface boundary. The unfilled circle represents the ghost point on an integer grid above the boundary while the filled circle indicates the mirror point on a fractional grid below the boundary. The filled squares represent the points on the integer grids below the surface while the unfilled squares indicate the points on the integer grids above the surface. The straight line across a pair of the mirror and ghost points is orthogonal to the free surface. In the immersed boundary method, the wavefield value of a ghost point is the negative value of the wavefield at the corresponding mirror point. 2.2. Computing the wavefield at mirror points iteratively The mirror points generally are at a fractional grid but the wavefield is computed at an integer grid using a conventional finite-difference method. To compute the wavefield value on the mirror points, we use three interpolation methods, namely bilinear interpolation, Lagrange interpolation and Kaiser windowed Sinc interpolation (Hicks 2002; Yao et al. 2018a; Li et al. 2019), rather than extrapolation because it could produce large errors, which is illustrated by an example in figure 2. Among the three interpolation methods, bilinear interpolation has the lowest accuracy because it ignores the curves of the waveform resulting in artificial high-frequency components, which is illustrated in figure 3. Figure 2. Open in new tabDownload slide Resampling a wavefield with Lagrange interpolation and extrapolation. The black dashed curve represents the original wavefield that is generated with a spatial interval of 1 m while the black cycles indicate the sampled points from the black dashed line with an interval of 10 m. The green, red and blue curves show the resampled wavefield back to the interval of 1 m by using symmetric interpolation, asymmetric interpolation and extrapolation respectively. The extrapolation uses four points on the left side of the target point while the symmetric interpolation takes four points from both the left and right side of the target point. The asymmetric interpolation takes four points on the left side and one point on the right side. Figure 2. Open in new tabDownload slide Resampling a wavefield with Lagrange interpolation and extrapolation. The black dashed curve represents the original wavefield that is generated with a spatial interval of 1 m while the black cycles indicate the sampled points from the black dashed line with an interval of 10 m. The green, red and blue curves show the resampled wavefield back to the interval of 1 m by using symmetric interpolation, asymmetric interpolation and extrapolation respectively. The extrapolation uses four points on the left side of the target point while the symmetric interpolation takes four points from both the left and right side of the target point. The asymmetric interpolation takes four points on the left side and one point on the right side. Figure 3. Open in new tabDownload slide Interpolation of one trace of a wavefield snapshot. The red dashed line indicates the original waveform with a sampling interval of 1 m. The sparse black circles represent the resample points at an interval of 10 m. The violet, blue and green curves indicate the results from the bilinear interpolation, Lagrange interpolation and Kaiser windowed Sinc interpolation, respectively. Figure 3. Open in new tabDownload slide Interpolation of one trace of a wavefield snapshot. The red dashed line indicates the original waveform with a sampling interval of 1 m. The sparse black circles represent the resample points at an interval of 10 m. The violet, blue and green curves indicate the results from the bilinear interpolation, Lagrange interpolation and Kaiser windowed Sinc interpolation, respectively. Figure 1 shows an example of interpolating the wavefield with 4 × 4 integer points. As can be seen, symmetric interpolation is used along both horizontal and vertical directions to compute the wavefield at each mirror point. However, the unknown ghost points also participate in the interpolation process, which causes the normal interpolation methods to become impractical. To solve this problem, we carry out the symmetric interpolation iteratively, which can be expressed as follows: $$\begin{equation} \left\{ {\begin{array}{@{}*{1}{c}@{}} {p_{mirror}^n{\rm{\, = }}\sum\limits_{i = 1}^m {{c_i}p_i^n\,} ,}\\ {p_{ghost}^{n + 1}\,\, = - p_{mirror}^n,} \end{array}} \right. \end{equation}$$(4) where |$p_i^n$| represents the pressure wavefield at the integer grid (including the ghost points participating the interpolation) in the n-th iteration, |${c_i}$| is the interpolation coefficients for the i-th point of interpolation, |$p_{mirror}^n$| denotes the pressure wavefield at the mirror points in the n-th iteration of interpolation and the second equation mirrors the wavefield of the mirror points to the corresponding ghost points antisymmetrically. In equation 4, |$p_i^n$| could include the pressure wavefield at ghost points. Our test shows no more than 20 iterations would converge the process of interpolation. 2.3. Calculating the particle vibration velocity at the surface In a land seismic survey, the geophone receivers are used to record the particle vibration velocity on the surface or very close to the surface as the pressure wavefield is zero at the surface. In addition, the surface is not flat, so that the position of the receivers is also on a fractional grid. To simulate the geophone record, we can use the following equation to convert the pressure wavefield near the surface into the particle vibration velocity: $$\begin{equation} \frac{{\partial {{\bf v}}({{\bf x}},t)}}{{\partial t}} = - \frac{1}{{\rho ({{\bf x}})}}\nabla p({{\bf x}},t), \end{equation}$$(5) where |${{\bf v}}$| represents the vector of particle vibration velocity, it has the horizontal and vertical components |${v_x}$| and |${v_z}$| for a 2D model, and |$\nabla $| indicates the operator of gradient. Afterward, we can use the interpolation methods mentioned above to compute the particle vibration velocity at the receiver locations. 2.4. Algorithm By combining all the theory and techniques described, the algorithm of our immersed boundary method for modelling irregular free surface is listed as follows: Setting the number of ghost layers to half the length of the finite-difference stencil; Calculating the coordinates of the intercept point for each ghost point by finding the closest point on the free surface to the ghost point; Computing the coordinates of mirror points; Setting the model parameters (e.g. velocity and density) of ghost points as those of the corresponding mirror points; Carrying out one time step of wavefield modelling by using conventional FDMs below the irregular surface; Setting the initial wavefield value of the ghost point to 0; Using a symmetric interpolation method to compute the wavefield value of the mirror points from the surrounding points at integer grids; Updating the wavefield value of ghost points with the negative value of mirror points; Returning to the step (7) to update the wavefield of mirror and ghost points iteratively. Usually, the update converges in 20 iterations; Calculating the particle vibration velocity in the region near the surface based on the obtained pressure wavefield by using equation 5, and interpolating the particle vibration velocity on the irregular surface; Returning to step (5) for the wavefield modelling of the next time step until the modelling finishes. 3. Numerical examples 3.1. A homogenous model with a Gaussian topography To verify the effectiveness of the proposed immersed boundary method, we used a Gaussian function to design an irregular surface model, which is shown in figure 4. The Gaussian function for the irregular surface is |$z(x) = 605 - 300{e^{ - {{( {\frac{{x - 1250}}{{200}}} )}^2}}}$|⁠, where the x represents the horizontal direction distance (units of meters) and |$z(x)$| indicates the depth of the surface. The model is discretised into 250 by 250 cells with a sampling interval of 10 m in both horizontal and depth directions. The source signature is a Ricker wavelet with a dominant frequency of 25 Hz. To generate the reference wavefield, we also discretised the model with a sampling interval of 1 meter (figure 4b), which is equivalent to 120 samples per wavelength of the dominant frequency, i.e. 25 Hz. Figure 4. Open in new tabDownload slide A homogenous model with a Gaussian topography. The discretised models with the intervals of (a) 10 m and (b) 1 m in both horizontal and vertical directions. Figure 4. Open in new tabDownload slide A homogenous model with a Gaussian topography. The discretised models with the intervals of (a) 10 m and (b) 1 m in both horizontal and vertical directions. Figure 5 shows the positional configuration of the ghost points and their mirror points in the irregular surface model. The black lines represent the irregular surface. The colored circles represent the four-layer ghost points while the colored filled circles indicate the mirror points corresponding to the ghost points. The small boxes at the bottom left and top right in figure 5 show the enlarged view. As can be seen, the ghost points are all on the integer grids while the mirror points are almost entirely on the fractional grids. Figure 5. Open in new tabDownload slide The positional configuration of the ghost points and their mirror points in the model with a Gaussian topography shown in figure 4. The black lines represent the irregular surface boundary and the unfilled color circles represent the four layers of ghost points, while the color-filled circles indicate the mirror points corresponding to every ghost point. In the two small boxes at the bottom left and top right are the enlarged views of the two sections indicated by the dashed line boxes. Figure 5. Open in new tabDownload slide The positional configuration of the ghost points and their mirror points in the model with a Gaussian topography shown in figure 4. The black lines represent the irregular surface boundary and the unfilled color circles represent the four layers of ghost points, while the color-filled circles indicate the mirror points corresponding to every ghost point. In the two small boxes at the bottom left and top right are the enlarged views of the two sections indicated by the dashed line boxes. The source is located in the center of the model along the horizontal direction and at a depth of 610 m. The receivers are spread evenly across the model at a depth of 610 m. In the test, we first used Lagrange interpolation for updating the wavefields of the ghost points. Figure 6 shows that the wavefield value at each ghost point on the first ghost layer against the number of iterations at the time of 300 ms. Each black line represents the evolution of the wavefield value at a ghost point. In the iteration process, the wavefield value at a ghost point is set to zero initially and then updated iteratively before reaching the preset maximum iteration number. As can be seen, the wavefield at the ghost points (the black stars) gradually converges to a stable state as increasing the number of iterations. Afterward, the next time step of the finite-difference modelling is carried out. Figure 6. Open in new tabDownload slide The wavefield value of the ghost points on the first ghost layer at the time of 300 ms against the number of iterations. Each black curve represents the wavefield at one position. Figure 6. Open in new tabDownload slide The wavefield value of the ghost points on the first ghost layer at the time of 300 ms against the number of iterations. Each black curve represents the wavefield at one position. For comparison, we used the vacuum method to simulate the wavefields with the spatial intervals of 10 m and 1 m in both directions. Four snapshots with the intervals of 10 m and 1 m are shown in figures 7a and b, respectively. As can be seen, a large number of diffraction waves are generated in the case of an interval of 10 m. The seismic records are shown in figures 8a and b. By comparison, the fine sampling interval can remove the artificial diffractions. However, it increases the computational cost significantly. For instance, the computational cost of the interval of 1 m is 100 times more than that of the interval of 10 m in this 2d case. Figure 7. Open in new tabDownload slide The seismic pressure wavefield snapshots using the vacuum method at the spatial intervals of: (a) 10 m and (b) 1 m, (c) our immersed boundary method with symmetric interpolation at a spatial interval of 10 m, (d) the immersed boundary method with asymmetric interpolation at a spatial interval of 10 m and (e) the immersed boundary method with extrapolation at a spatial interval of 10 m. The recording times from left to right are 150, 300, 450 and 600 ms. Note that the ghost wavefield is not shown. Figure 7. Open in new tabDownload slide The seismic pressure wavefield snapshots using the vacuum method at the spatial intervals of: (a) 10 m and (b) 1 m, (c) our immersed boundary method with symmetric interpolation at a spatial interval of 10 m, (d) the immersed boundary method with asymmetric interpolation at a spatial interval of 10 m and (e) the immersed boundary method with extrapolation at a spatial interval of 10 m. The recording times from left to right are 150, 300, 450 and 600 ms. Note that the ghost wavefield is not shown. Figure 8. Open in new tabDownload slide The profiles of seismic records using the vacuum method at the spatial intervals of: (a) 10 m and (b) 1 m, (c) our immersed boundary method with symmetric interpolation at a spatial interval of 10 m and (d) the immersed boundary method with asymmetric interpolation at a spatial interval of 10 m. Figure 8. Open in new tabDownload slide The profiles of seismic records using the vacuum method at the spatial intervals of: (a) 10 m and (b) 1 m, (c) our immersed boundary method with symmetric interpolation at a spatial interval of 10 m and (d) the immersed boundary method with asymmetric interpolation at a spatial interval of 10 m. As can be seen, the vacuum method with a fine grid can remove artificial diffractions presented in figure 7a at the cost of significant computational increase. By contrast, our immersed boundary method with an interval of 10 m also achieves a high-quality result shown in figure 7c. This improvement also can be seen in the record profiles shown in figure 8c. For comparison, we used asymmetric interpolation to calculate the forward wavefield. The configuration is the same as figure 4 of Hu (2016). We used iterative interpolation (Equation 7 of Hu 2016), which is similar to Equation 4, to compute the wavefield at fictitious points denoted with triangles in figure 4 of Hu (2016). Then, the wavefield value at one mirror point is computed by using asymmetric interpolation with the wavefield values at four fictitious points and one point on the free surface, the wavefield value of which is zero. As shown in figure 7d, the modelling becomes unstable after 300 ms. The corresponding shot record is shown in figure 8d. In addition, we also tried to compute the wavefield at ghost points by directly using extrapolation with the wavefield value at four fictitious points and one point on the free the free surface, which is the method illustrated by equations 4–7 by Hu (2016). The extrapolation coefficients in this example are computed using Lagrange polynomials. This way leads the modelling to unstable even faster, which is shown in figure 7e. Extrapolation generating large errors is also demonstrated in figure 2. This example demonstrates that the immersed boundary method requires high accuracy of calculation for wavefields at the mirror points or ghost points. Our iterative symmetric interpolation is a choice for this purpose. To further illustrate the convergence of the iterative interpolation method, we plotted out a single trace at a distance of 1000 m. Figure 9 parts a and b show the seismic trace record received with different iterations numbers during computing wavefields of ghost points. As can be seen, the artificial diffractions have been reduced after the first iteration. Besides, the wavefield converges to a stable state quickly by increasing the iterations number. With 20 iterations, the result is very similar to the record generated with a fine sampling interval (dx = dz = 1 m). Figure 9. Open in new tabDownload slide The seismograms of the trace at a distance of 1000 m. (a, b) The seismograms with different iterations. (c, d) The seismograms with different interpolation methods at an iteration of 20. The dashed curves represent the result using vacuum methods. (b) and (d) show the enlarged view of a section of the trace indicated by a dashed box in (a) and (c), respectively. Figure 9. Open in new tabDownload slide The seismograms of the trace at a distance of 1000 m. (a, b) The seismograms with different iterations. (c, d) The seismograms with different interpolation methods at an iteration of 20. The dashed curves represent the result using vacuum methods. (b) and (d) show the enlarged view of a section of the trace indicated by a dashed box in (a) and (c), respectively. To verify the influence of different interpolation algorithms for our immersed boundary method, we test the three interpolation methods, i.e. bilinear interpolation, Lagrange interpolation and Kaiser windowed Sinc interpolation. The seismogram of the result is shown in figure 9 parts c and d. By comparison, it is clear that Lagrange interpolation and Kaiser windowed Sinc interpolation give similar results that were much better than bilinear interpolation. Since Lagrange interpolation has a lower computational cost than Kaiser windowed Sinc interpolation, we show the results of the rest tests using it. In a real land seismic survey, the receivers are usually plugged on the surface instead of underground to record the particle vibration velocity. In our modelling method, we generate the surface record by using equation 5 together with the Kaiser windowed Sinc interpolation. The receivers are placed on the nearest integer points for the vacuum method. The surface record comparison between the vacuum method and our method was shown in figure 10, which corresponds to the wavefield snapshots shown in figure 7a–c. It is obvious that the results from the vacuum method have errors, which are proportional to the grid sampling interval. The larger the grid interval, the worse the results from the vacuum method. By contrast, our method shows better results. Figure 10. Open in new tabDownload slide The seismic profiles of particle vibration velocity on the surface. The vertical components are shown in the left column while the horizontal components are in the right column. Three rows from top to bottom are generated using the vacuum method with an interval of 10 and 1 m and our immersed boundary method at an interval of 10 m, which corresponds to the wavefield snapshots in parts (a–c), respectively. Figure 10. Open in new tabDownload slide The seismic profiles of particle vibration velocity on the surface. The vertical components are shown in the left column while the horizontal components are in the right column. Three rows from top to bottom are generated using the vacuum method with an interval of 10 and 1 m and our immersed boundary method at an interval of 10 m, which corresponds to the wavefield snapshots in parts (a–c), respectively. We also put the source at 5 m below the peak of the Gaussian hill and record the particle vibration velocity on the surface. The results are shown in figure 11, which shows clear direct arrivals. Figure 11. Open in new tabDownload slide The seismic profiles of the particle vibration velocity: (a) the horizontal component and (b) the vertical component. The source is located at 5 m below the peak of the Gaussian hill while the receivers are at the surface. Figure 11. Open in new tabDownload slide The seismic profiles of the particle vibration velocity: (a) the horizontal component and (b) the vertical component. The source is located at 5 m below the peak of the Gaussian hill while the receivers are at the surface. 3.2. The homogenous models with different sine-shaped topography In the second example, we used four sine functions to design irregular surface models shown in figure 12. The size of the model and discretisation parameters are the same for the previous model. The sine function for defining the topography is |$z(x) = 405 - 100\sin ( {\frac{{\pi x}}{a}} )$|⁠, where a determines the frequency of topography oscillation. a is set to 100, 150, 200 and 300 in these four models shown in figure 12 parts a–d, respectively. Consequently, the maximum dip angles for the four models are about 73, 64, 57 and 46°. The source wavelet is a 30 Hz Ricker wavelet. Figure 12. Open in new tabDownload slide Four homogeneous models with different sine-shaped topography, which represent different steepnesses and complexities of the surface. The topography of the four models is defined by the sine functions |$z(x) = 405 - 100\sin ( {\frac{{\pi x}}{a}} )$|⁠, where the values of a are 100, 150, 200 and 300 in (a), (b), (c) and (d), respectively. These models are sampled with a spatial interval of 10 m in both horizontal and vertical directions. Figure 12. Open in new tabDownload slide Four homogeneous models with different sine-shaped topography, which represent different steepnesses and complexities of the surface. The topography of the four models is defined by the sine functions |$z(x) = 405 - 100\sin ( {\frac{{\pi x}}{a}} )$|⁠, where the values of a are 100, 150, 200 and 300 in (a), (b), (c) and (d), respectively. These models are sampled with a spatial interval of 10 m in both horizontal and vertical directions. Figure 13 shows the records generated by setting the source at the center along the horizontal direction and a depth of 610 m and the receivers at a depth of 610 m for these four models. The four rows from top to bottom in figure 13 depict the records from the four models shown in figure 12 parts a–d, respectively. The left and middle columns of figure 13 show the records by using the vacuum method with the spatial sampling intervals of 10 and 1 m, respectively. The right column shows the records by using our immersed boundary method at spatial sampling intervals of 10 m in both horizontal and vertical directions. As shown in figure 13, the vacuum method with a sampling interval of 10 m generates artificial diffractions for all four models. By contrast, our immersed boundary method can eliminate these diffractions effectively even for the model shown in figure 12a, which has a maximum dip angle of 73° and strongly bending topography. Note that our immerse boundary method is stable for all the four models (we tested modelling for a duration of 10 s). Figure 13. Open in new tabDownload slide The recorded seismic profiles from the four sine-shaped topography models shown in figure 12. The four rows from top to bottom depict the records from the four models shown in (a–d), respectively. The left and middle columns show the records by using the vacuum method with the spatial sampling intervals of 10 and 1 m, respectively. The right column shows the records by using our immersed boundary method at a spatial sampling interval of 10 m in both horizontal and vertical directions. The red arrows highlight the artificial scattering in wavefields for each group of records. Figure 13. Open in new tabDownload slide The recorded seismic profiles from the four sine-shaped topography models shown in figure 12. The four rows from top to bottom depict the records from the four models shown in (a–d), respectively. The left and middle columns show the records by using the vacuum method with the spatial sampling intervals of 10 and 1 m, respectively. The right column shows the records by using our immersed boundary method at a spatial sampling interval of 10 m in both horizontal and vertical directions. The red arrows highlight the artificial scattering in wavefields for each group of records. If we reduce the value of a more, which gives steeper dip angles and sharper topography, our immersed boundary method starts to generate noticeable diffractions and become even unstable. This is because that very sharp topography causes many mirror points in one cell and the wavefield energy in one cell is spread to several ghost points. In fact, it is almost impractical to carry out a seismic survey in an area where the surface slope is very large due to the difficulty of moving the equipment. As a result, our method should be good enough to handle the real topography problems for seismic exploration. 3.3. The Canadian foothill model Finally, we applied our immersed boundary method to a more realistic model, the Canadian foothill model, which has complicated surface topography (figure 14). The model has a size of 1668 cells by 1000 cells in the horizontal and vertical directions with a sampling interval of 15 and 10 m, respectively. A 25 Hz Ricker wavelet is used as the source signature. The time sampling interval is 0.5 ms. One shot record in the middle of the model is generated. The source is located on the first integer cell below the surface indicated by the yellow star in figure 14. The receivers are deployed on the surface indicated by the dashed red line in figure 14. Figure 14. Open in new tabDownload slide The velocity profile of the Canadian foothill model. The model is discretised into 1668 cells by 1000 cells with the spatial intervals of 15 m and 10 m in the horizontal and vertical directions, respectively. The dashed red line indicates the surface of the model. The yellow star represents the source location of the test shot. The receiver array is located on the surface. Figure 14. Open in new tabDownload slide The velocity profile of the Canadian foothill model. The model is discretised into 1668 cells by 1000 cells with the spatial intervals of 15 m and 10 m in the horizontal and vertical directions, respectively. The dashed red line indicates the surface of the model. The yellow star represents the source location of the test shot. The receiver array is located on the surface. The snapshots of the pressure wavefield at 1000, 1500 and 2500 ms are shown in figure 15, respectively. We can see the reflected waves from the free surface. The particle vibration velocity profiles of the shot are shown in figure 16. It shows that our immersed boundary method is stable. As can be seen from figure 16, the seismic energy distribution is uneven due to the influence of the irregular topography. In the profile of the horizontal component, the events’ amplitudes are reduced to near zero at the position where the surface is close to horizontal (figure 16a). However, these events are clear and continuous in the profile of the vertical component (figure 16b). This is the main reason why only the vertical component is acquired in a land seismic data survey usually. It can also be seen that the amplitudes of events are attenuated if a valley sits between the source and receivers. Figure 15. Open in new tabDownload slide The snapshots of the pressure wavefield at (a) 1000 ms, (b) 1500 ms and (c) 2500 ms. The velocity is plotted as the background. Figure 15. Open in new tabDownload slide The snapshots of the pressure wavefield at (a) 1000 ms, (b) 1500 ms and (c) 2500 ms. The velocity is plotted as the background. Figure 16. Open in new tabDownload slide The seismic profiles of the horizontal (a) and vertical (b) components of particle vibration velocity. Figure 16. Open in new tabDownload slide The seismic profiles of the horizontal (a) and vertical (b) components of particle vibration velocity. By calculating the particle vibration velocity on the surface at fractional grids, this example demonstrates our method could be used in actual land seismic wavefield modelling. The numerical example on the Canadian foothill model, which contains a complex topography, also validates that our immersed boundary method with iterative symmetric interpolation has applicability for a realistic land seismic data modelling. 4. Conclusion In this work, we proposed an immersed boundary method with iterative symmetric interpolation. This method can effectively handle the irregular topography for acoustic seismic wavefield modelling on rectangular meshes. As the rectangular meshes are used for the finite-difference modelling, the merits of stability and efficiency of FDMs remain. In this method, first, the number of ghost layers is set to half the length of the finite-difference stencil. Second, the mirror points are found by orthogonally projecting the ghost points to fractional points below the free surface. Third, the wavefield at these mirror points is calculated by using an iterative symmetric interpolation method. Finally, the wavefield at these ghost points is the negative value of the wavefield at the corresponding mirror points. Since the number of ghost layers can be set freely, this method can be used for any order finite-difference modelling. The iterative symmetric interpolation can allow the ghost points to participate in the interpolation; consequently, the symmetric interpolation with uniformly sampled points is used to ensure the accuracy of the ghost wavefield and furthermore the stability of finite-difference modelling. Our numerical tests show that both Lagrange and Kaiser windowed Sinc interpolation can give a high-quality ghost wavefield. To simulate the land seismic acquisition, we convert the pressure wavefield into the particle vibration velocity at the free surface. The numerical tests on three models demonstrate the effectiveness of this immersed boundary method. Therefore, it is ready for seismic imaging and inversion with irregular surface topography. Acknowledgements This work was partially supported by National Key R&D Program of China (grant no. 2018YFA0702502), the National Natural Science Foundation of China (grant no. 41974142), National Science Foundation (grant no. EAR-1547228) and Science Foundation of China University of Petroleum (Beijing) (grant no. 2462019YJRC005). Conflict of interest statement. None declared. References Almuhaidib A.M. , Toksoz M.N., 2015 . Finite difference elastic wave modelling with an irregular free surface using ADER scheme , Journal of Geophysics and Engineering , 12 , 435 – 447 . Google Scholar Crossref Search ADS WorldCat Alterman Z.S. , Karal F.C.J., 1968 . Propagation of elastic waves in layered media by finite difference methods , Bulletin of the Seismological Society of America , 58 , 367 – 398 . OpenURL Placeholder Text WorldCat Beasley C.J. , Lynn W., 1989 , Zero-velocity layer: migration from irregular surfaces , Geophysics , 22 , 35 – 40 . OpenURL Placeholder Text WorldCat Bleibinhaus F. , Rondenay S., 2009 . Effects of surface scattering in full-waveform inversion , Geophysics , 74 , 69 – 77 . Google Scholar Crossref Search ADS WorldCat Bohlen T. , Saenger E.H., 2006 . Accuracy of heterogeneous staggered-grid finite-difference modelling of Rayleigh waves , Geophysics , 71 , T109 – T115 . Google Scholar Crossref Search ADS WorldCat Borisov D. , Modrak R., Gao F.C., Tromp J., 2018 . 3D elastic full-waveform inversion of surface waves in the presence of irregular topography using an envelope-based misfit function , Geophysics , 83 , R1 – R11 . Google Scholar Crossref Search ADS WorldCat Cao J. , Chen J.B., 2018 . A parameter-modified method for implementing surface topography in elastic-wave finite-difference modelling , Geophysics , 83 , T313 – T332 . Google Scholar Crossref Search ADS WorldCat Gao L.F. , Brossier R., Pajot B., Tago J., Virieux J., 2015 . An immersed free-surface boundary treatment for seismic wave simulation , Geophysics , 80 , T193 – T209 . Google Scholar Crossref Search ADS WorldCat Graves R.W. , 1996 . Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences , Bulletin of the Seismological Society of America , 86 , 1091 – 1106 . OpenURL Placeholder Text WorldCat Hayashi K. , Burns D.R., Toksoz M.N., 2001 . Discontinuous-grid finite-difference seismic modelling including surface topography , Bulletin of the Seismological Society of America , 91 , 1750 – 1764 . Google Scholar Crossref Search ADS WorldCat He Y. , Wang H.Z., Ma Z.T., Hu Z.B, 2002 . Pre-stack wave equation depth migration for irregular topography , Progress Exploration Geophysics -CH , 25 , 13 – 19 . OpenURL Placeholder Text WorldCat Hicks G.J. , 2002 . Arbitrary source and receiver positioning in finite-difference schemes using Kaiser windowed sinc functions , Geophysics , 67 , 156 – 166 . Google Scholar Crossref Search ADS WorldCat Hu W.Y. , 2016 . An improved immersed boundary finite-difference method for seismic wave propagation modelling with arbitrary surface topography , Geophysics , 81 , T311 – T322 . Google Scholar Crossref Search ADS WorldCat Jastram C. , Tessmer E., 2010 . Elastic modelling on a grid with vertically varying spacing , Geophysical Prospecting , 42 , 357 – 370 . Google Scholar Crossref Search ADS WorldCat Kelly K.R. , Ward R.W., Treitel S., Alford R.M., 1976 . Synthetic seismograms: a finite difference approach , Geophysics , 41 , 2 – 27 . Google Scholar Crossref Search ADS WorldCat Komatitsch D. , Vilotte J.P., 1998 . The spectral-element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures , Bulletin of the Seismological Society of America , 88 , 368 – 392 . OpenURL Placeholder Text WorldCat Kosloff D.D. , Baysal E., 1982 . Forward modelling by a Fourier method , Geophysics , 47 , 1402 – 1412 . Google Scholar Crossref Search ADS WorldCat Levander A.R. , 1988 . Fourth-order finite-difference P-SV seismograms , Geophysics , 53 , 1425 – 1436 . Google Scholar Crossref Search ADS WorldCat Li Q.Y. , Wu G.C., Wu J.L., Duan P., 2019 . Finite difference seismic forward modelling method for fluid-solid coupled media with irregular seabed interface , Journal of Geophysics and Engineering , 16 , 198 – 214 . Google Scholar Crossref Search ADS WorldCat Li Y.Y. , Li Z.C., Zhang K., Lin Y.Z., 2016 . Frequency-domain full waveform inversion with rugged free surface based on variable grid finite-difference method , Journal of Seismic Exploration , 25 , 543 – 559 . OpenURL Placeholder Text WorldCat Liu S.L. , Li X.F., Wang W.S., Liu Y.S., 2014 . A mixed-grid finite element method with PML absorbing boundary conditions for seismic wave modelling , Journal of Geophysics and Engineering , 11 , 1 – 13 . OpenURL Placeholder Text WorldCat Lombard B. , Piraux J., 2004 . Numerical treatment of two-dimensional interfaces for acoustic and elastic waves , Journal of Computational Physics , 195 , 90 – 116 . Google Scholar Crossref Search ADS WorldCat Lombard B. , Piraux J., Gélis C., Virieux J., 2008 . Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves , Geophysical Journal International , 172 , 252 – 261 . Google Scholar Crossref Search ADS WorldCat Marfurt K.J. , 1984 . Accuracy of finite-difference and finite-element modelling of the scalar and elastic wave equations , Geophysics , 49 , 533 – 549 . Google Scholar Crossref Search ADS WorldCat Meng W.J. , Fu L.Y., 2017 . Seismic wavefield simulation by a modified finite element method with a perfectly matched layer absorbing boundary , Journal of Geophysics and Engineering , 14 , 852 – 864 . Google Scholar Crossref Search ADS WorldCat Mittal R. , Dong H., Bozkurttas M., Najjar F.M., Vargas A., von Loebbecke A., 2008 . A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries , Journal of Computational Physics , 227 , 4825 – 4852 . Google Scholar Crossref Search ADS PubMed WorldCat Muir F. , Dellinger J., Etgen J., Nichols D., 1992 . Modelling elastic fields across irregular boundaries , Geophysics , 57 , 1189 – 1193 . Google Scholar Crossref Search ADS WorldCat Mulder W.A. , Huiskes M.J., 2017 . A simple finite-difference scheme for handling topography with the first-order wave equation , Geophysical Journal International , 210 , 482 – 499 . Google Scholar Crossref Search ADS WorldCat Nielsen P. , Flemming L., Berg P., Skovgaard O., 1994 . Using the pseudospectral technique on curved grids for 2D acoustic forward modelling , Geophysical Prospecting , 42 , 321 – 341 . Google Scholar Crossref Search ADS WorldCat Nuber A. , Manukyan E., Maurer H., 2016 . Ground topography effects on near-surface elastic full waveform inversion , Geophysical Journal International , 207 , 67 – 71 . Google Scholar Crossref Search ADS WorldCat Oprsal I. , Zahradník J., 1999 . Elastic finite‐difference method for irregular grids , Geophysics , 64 , 240 – 250 . Google Scholar Crossref Search ADS WorldCat Oral E. , Gelis C., Bonilla L.F., 2019 . 2-D P-SV and SH spectral element modelling of seismic wave propagation in non-linear media with pore-pressure effects , Geophysical Journal International , 217 , 1353 – 1365 . Google Scholar Crossref Search ADS WorldCat Peskin C.S. , 1972 . Flow patterns around heart valves: a numerical method , Journal of Computational Physics , 10 , 252 – 271 . Google Scholar Crossref Search ADS WorldCat Pratt R. , 1999 . Seismic waveform inversion in the frequency domain, Part 1: theory and verification in a physical scale model , Geophysics , 64 , 888 – 901 . Google Scholar Crossref Search ADS WorldCat Rao Y. , Wang Y., 2013 . Seismic waveform simulation with pseudo-orthogonal grids for irregular topographic models , Geophysical Journal International , 194 , 1778 – 1788 . Google Scholar Crossref Search ADS WorldCat Rao Y. , Wang Y., 2018 . Seismic waveform simulation for models with fluctuating interfaces , Scientific Reports , 8 , 3098 . Google Scholar Crossref Search ADS PubMed WorldCat Reshef M. , 1991 . Depth migration from irregular surfaces with depth extrapolation methods , Geophysics , 56 , 119 – 122 . Google Scholar Crossref Search ADS WorldCat Shtivelman V. , Canning A., 1988 . Datum correction by wave equation extrapolation , Geophysics , 53 , 1311 – 1322 . Google Scholar Crossref Search ADS WorldCat Sun Y.C. , Zhang W., Chen X., 2016 . Seismic-wave modelling in the presence of surface topography in 2D general anisotropic media by a curvilinear grid finite-difference method , Bulletin of the Seismological Society of America , 106 , 1036 – 1054 . Google Scholar Crossref Search ADS WorldCat Sun Y.C. , Zhang W., Xu J.K., Chen X., 2017 . Numerical simulation of 2-D seismic wave propagation in the presence of a topographic fluid-solid interface at the sea bottom by the curvilinear grid finite-difference method , Geophysical Journal International , 210 , 1721 – 1738 . Google Scholar Crossref Search ADS WorldCat Tarantola A. , 1984 . Inversion of seismic reflection data in the acoustic approximation , Geophysics , 49 , 1259 – 1266 . Google Scholar Crossref Search ADS WorldCat Tessmer E. , Kosloff D., Behle A., 1992 . Elastic wave propagation simulation in the presence of surface topography , Geophysical Journal International , 108 , 621 – 632 . Google Scholar Crossref Search ADS WorldCat Tseng Y.H. , Ferziger J.H., 2003 . A ghost-cell immersed boundary method for flow in complex geometry , Journal of Computational Physics , 192 , 593 – 623 . Google Scholar Crossref Search ADS WorldCat Virieux J. , 1986 . P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method , Geophysics , 51 , 1933 – 1942 . Google Scholar Crossref Search ADS WorldCat Wang Z.Y. , Huang J.P., Liu D.J., Li Z.C., Yong P., Yang Z.J., 2019 . 3D variable-grid full-waveform inversion on GPU , Petroleum Science , 16 , 1001 – 1014 . Google Scholar Crossref Search ADS WorldCat Warner M. et al. ., 2013 . Anisotropic 3D full-waveform inversion , Geophysics 78 , R59 – R80 . Google Scholar Crossref Search ADS WorldCat Yao G. , da Silva N.V., Debens H.A., Wu D., 2018 . Accurate seabed modelling using finite difference methods , Computational Geosciences , 22 , 469 – 484 . Google Scholar Crossref Search ADS WorldCat Yao G. , da Silva N.V., Wu D., 2018 . An effective absorbing layer for the boundary condition in acoustic seismic wave simulation , Journal ofGeophysics and Engineering , 15 , 495 – 511 . Google Scholar Crossref Search ADS WorldCat Yao G. , da Silva N.V., Wu D., 2019 . Reflection-waveform inversion regularized with structure-oriented smoothing shaping , Pure and Applied Geophysics , 176 , 5315 – 5335 . Google Scholar Crossref Search ADS WorldCat Yao G. , Wu D., Debens H.A., 2016 . Adaptive finite difference for seismic wavefield modelling in acoustic media , Scientific Reports , 6 , 30302 . Google Scholar Crossref Search ADS PubMed WorldCat Zahradník J. , Moczo P., Hron H., 1993 . Testing four elastic finite-difference schemes for behavior at discontinuities , Bulletin of the Seismological Society of America , 83 , 107 – 129 . OpenURL Placeholder Text WorldCat Zhang W. , Chen X.F., 2006 . Traction image method for irregular free surface boundaries in finite difference seismic wave simulation , Geophysical Journal International , 167 , 337 – 353 . Google Scholar Crossref Search ADS WorldCat Zhang W. , Shen Y., Zhao L., 2012 . Three-dimensional anisotropic seismic wave modelling in spherical coordinates by a collocated-grid finite-difference method , Geophysical Journal International , 188 , 1359 – 1381 . Google Scholar Crossref Search ADS WorldCat Zhao S. , 2010 . A fourth order finite difference method for waveguides with curved perfectly conducting boundaries , Computer Methods in Applied Mechanics and Engineering , 199 , 2655 – 2662 . Google Scholar Crossref Search ADS WorldCat © The Author(s) 2020. Published by Oxford University Press on behalf of the Sinopec Geophysical Research Institute. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
TI - An immersed boundary method with iterative symmetric interpolation for irregular surface topography in seismic wavefield modelling
JF - Journal of Geophysics and Engineering
DO - 10.1093/jge/gxaa019
DA - 2020-08-01
UR - https://www.deepdyve.com/lp/oxford-university-press/an-immersed-boundary-method-with-iterative-symmetric-interpolation-for-LZ1IGMmxDZ
SP - 643
EP - 660
VL - 17
IS - 4
DP - DeepDyve
ER -