An immersed boundary method with iterative symmetric interpolation for irregular surface topography in seismic wavefield modelling

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.


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. 2016Sun et al. , 2017 and the acoustic wavefield (Rao & Wang 2013 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. 644 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.

Methodology
The acoustic wave equation can be expressed as where p(x, t) represents the pressure field at the location x and the time t, s(t) denotes the source signature, v(x) and (x) are the velocity and density respectively, (x − x s ) represents the Dirichlet delta function for injecting the point source at the source location 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, and the free surface boundary condition, where ∑ 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).

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.

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 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: where p n i 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 n mirror 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 n i could include the pressure wavefield at ghost points. Our test shows no more than 20 iterations would converge the process of interpolation.

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: where 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 ∇ indicates the operator of gradient.  Afterward, we can use the interpolation methods mentioned above to compute the particle vibration velocity at the receiver locations.

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: (1) Setting the number of ghost layers to half the length of the finite-difference stencil; (4) Setting the model parameters (e.g. velocity and density) of ghost points as those of the corresponding mirror points; (5) Carrying out one time step of wavefield modelling by using conventional FDMs below the irregular surface; (6) Setting the initial wavefield value of the ghost point to 0; (7) Using a symmetric interpolation method to compute the wavefield value of the mirror points from the surrounding points at integer grids; (8) Updating the wavefield value of ghost points with the negative value of mirror points; (9) Returning to the step (7) to update the wavefield of mirror and ghost points iteratively. Usually, the update converges in 20 iterations; field by using equation 5, and interpolating the particle vibration velocity on the irregular surface; (11) Returning to step (5) for the wavefield modelling of the next time step until the modelling finishes.

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( 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 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.
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 finitedifference modelling is carried out.
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.
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 demon-strated 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  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). 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.
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.

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( 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 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).
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.   Figure 13. 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.

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

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.