Semi-exact local absorbing boundary condition for seismic wave simulation

An absorbing boundary condition is necessary in seismic wave simulation for eliminating the unwanted artificial reflections frommodel boundaries. Existing boundary condition methods often have a trade-off between numerical accuracy and computational efficiency. We proposed a local absorbing boundary condition for frequency-domain finite-difference modelling. The proposed method benefits from exact local plane-wave solution of the acoustic wave equation along predefined directions that effectively reduces the dispersion in other directions. This method has three features: simplicity, accuracy and efficiency. Numerical simulation demonstrated that the proposed method has higher efficiency than the conventional methods such as the second-order absorbing boundary condition and the perfectly matched layer (PML) method. Meanwhile, the proposed method shared the same low-cost feature as the first-order absorbing boundary condition method.


Introduction
Finite-difference solution of acoustic wave equation is an engine of advanced imaging technologies such as full waveform inversion (FWI) and reverse time migration. Frequencydomain finite-difference (FDFD) methods are of special interest, as they offer advantages over the time-domain counterparts. The advantages include simultaneous multisource implementation (Wang 2011), easier modelling of absorptive media Amini & Javaherian 2011) and selective control over the frequencies ( Jo et al. 1996;).
Because of limited computer memory available, an unbounded domain needs to be trimmed for numerical seismic wave simulation. As a result, it is imperative to design an absorbing boundary condition to absorb incoming waves toward the numerical boundaries (Zhao et al. 2019). Three commonly used absorbing boundary conditions in seismic wave simulation are: one-way wave equation absorbing boundary condition (ABC), absorbing sponge layer (ASL) and perfectly matched layer (PML).
The conventional ABC method is a low-order approximation of a one-way wave equation (Clayton & Engquist 1977;Higdon 1986). A two-way wave equation may be decomposed into two one-way wave equations in two opposite directions, called upgoing and downgoing directions. Because the ABC method is a local boundary condition, it is computationally the most efficient absorbing boundary method. However, the ABC shows good performance only when the outgoing wave propagation angle is nearly normal to the boundary. Addressing this drawback requires employing high-order approximations (Xu & McMechan 2012;Ren & Liu 2013;Liu & Sen 2018). Then, a high-order ABC method requires high-order spatial derivatives, which are inconvenient in practical implementations (Givoli 2004;Rabinovich et al. 2010;Song et al. 2015).
The ASL method attempts to eliminate the incident waves by introducing attenuation inside a padded layer. The attenuation is realised through a gradually boosted damping function inside the padded layer (Cerjan et al. 1985;Shin 1995;Hall & Wang 2012). Two challenges of the ASL method are finding an optimum width of the padded layer and a proper damping function. The ASL method can be more effective than a low-order ABC method, if the damping layer is relatively thick with optimal attenuation parameters in the damping function (Sochacki et al. 1987;Zhou 1988). However, too thick a padded layer will lead to extra cost in running time.
The PML method, like the ASL method, introduces the attenuation inside a padded layer to the model. However, a PML method incorporates the attenuation via a physically consistent approach (Berenger 1994;Collino & Tsogka 2001;Hustedt et al. 2004;Rao et al. 2016;He et al. 2020). The basic idea is to modify the partial derivatives in the wave equation through coordinate stretching (Abarbanel & Gottlieb 1997). In a continuous domain, the PML method is proved to be effective with completely free of boundary reflections. However, a finite length of the padded layer and the discretisation downgrade its effectiveness and result in angle and frequency dependencies (Collino & Tsogka 2001). The phenomena are manifested as unwanted reflections from boundaries particularly for low-frequency waves, the evanescent waves and the nearly grazing incident waves (Gao et al. 2017). An unsplit convolutional PML method could tackle these drawbacks in the classic PML method (Komatitsch & Martin 2007;Martin & Komatitsch 2009).
The PML method has been widely used in seismic wave simulation, but the computational cost becomes prohibitive in the large-scale 3D applications, in which cases the cost is increased exponentially as the number of grids increases (Amini & Javaherian 2011). In this paper, we develop an efficient and low-cost ABC for large-scale modelling applications. Our development is inspired by the local plane-wave solution of the wave equation, which is an accurate and cost-effective FDFD scheme for solving the acoustic wave equation (Izadian et al. 2019). The basic idea of the local plane-wave solution is to force the dispersion to be zero in number of predefined directions that results in a more accurate finite-difference (FD) scheme. Our contribution here is to explore this idea for the ABC. We use a plane-wave solution of the acoustic wave equation to interpolate its dispersion relation in half of the directions available (between ± ∕2) and leave the rest of directions to be zero. This automatically gives a dispersion relation equivalent to that of the one-way wave equation.
In this paper, we first review the theory of the plane-wave interpolation for acoustic wave modelling, and then design a local ABC. We examine the performance of the proposed absorbing boundary condition against representatives of the ABC and PML method via various numerical examples.

The dispersion relation and approximations
To design an absorbing boundary condition based on planewave interpolation, it is inspirational to look at the ABC method. The 2D acoustic wave equation in the frequency domain is 2 P (x, z, ) where P(x, z, ) is the pressure wavefield, and is the angular frequency. Taking the spatial Fourier transform of equation (1) yields the following dispersion relation: where k x and k z are the horizontal and vertical wavenumbers, respectively. Equation (2) may be decomposed into two parts (Clayton & Engquist 1977): where the positive sign denotes the downgoing wave (positive z direction points downward), and the negative sign denotes the upgoing wave. Thus, for the top boundary, we have the following equation: The ABC method for the top boundary is an approximation to equation (4). The first-and second-order Taylor expansions are (Clayton & Engquist 1977) vk z = 1, The third-order Padé approximation is (Clayton & Engquist 1977) As discussed, the decomposition of the dispersion relation of acoustic wave equation has removed a half of the wavefield and, then, the ABC method designs the absorbing boundary condition by approximating the remaining half of wavefield.
In this paper, we propose a fundamentally different way to directly approximate the dispersion relation of the one-way wave equation (e.g. only the upgoing direction) by a planewave interpolation (PWI) solution. In the next section, we describe interpolation of the full dispersion relation which gives a highly accurate FD scheme to solve the acoustic wave equation. After that, we extend the concept to design a local ABC.

Wave modelling using plane-wave interpolation
The wave equation (1) may be rewritten in a general form of the Laplacian operator as ( Jo et al. 1996): where ∇ 2 0 and ∇ 2 45 denote the Cartesian and 45°rotated Laplacian operators, respectively (figure 1). Assuming P m,n = P(mh, nh) and h is the grid interval and considering the general form of the FD approximations of the Cartesian and rotated Laplacian operators: Plugging equations (9) and (10) into equation (8) gives with For any arbitrary direction, the plane-wave solution is P = e ik(x sin +zcos ) , where i is the imaginary symbol, and is the propagation angle with respect to the z-axis. Substituting the solution into equation (11) gives For the nine-point scheme, equation (13) is used for waves propagating along the Cartesian coordinate directions ( = 0, ± 2 , ) as well as a 45 0 rotated direction ( = ± 4 , ± 3 4 ). In this way, dispersion errors along these directions becomes zero.
where r = cos( 4 ) . From these two equations, we obtain coefficients C 1 and C 2 as follows: Finally, we substitute the coefficients into equation (13) and obtain a new nine-point scheme: This is an improved version of the standard nine-point scheme with more accuracy. The improvement is achieved by taking advantage of smoothness of the dispersion function and zeroing out dispersion along predefined directions, which consequently reduces error along other directions (Izadian et al. 2019). Figure 2 compares dispersion curves of the optimised nine-point ( Jo et al. 1996), optimised 25-point  (Shin & Sohn 1998) and PWI schemes for different propagation angles ( = 0, 12 , 6 , 4 ). It is evident that the PWI method outperforms than 25-point scheme in terms of numerical dispersion. PWI has constrained the dispersion error to be less than 0.4% for 2.5 grids per minimum wavelength.

Local absorbing boundary condition
We insert plane waves with a limited range of directions into a proper general form of a FD scheme for the ABC. This process will automatically interpolate only a limited range of directions and will keep the rest of the directions to zero. Therefore, the resulting FD scheme will propagate the waves only in that range of directions. In the following, we describe the ABCs for the top boundary and the top-left corner. These boundary conditions may be extended straightforwardly to other edges and corners.
Applying equation (18) along = (0, 4 , 2 ), we obtain corresponding equations for this stencil: c 1 + c 2 + e ikh c 3 + e ikh c 4 = 0, c 1 + e ikh c 2 + c 3 + e ikh c 4 = 0, c 1 + e jkrh c 2 + e jkrh c 3 + e 2jkrh c 4 = 0, where r = cos( 4 ) . Solving this system of equations in the same way as solving equation (20), we arrive to the following scheme: P m,n + P m+1,n + P m,n+1 + P m+1,n+1 = 0, Note that equations (22) and (25) came out of a fundamentally different approach as an ABC method. We started with a FD scheme and forced it to follow the exact twoway dispersion relation in some directions. This method exploits the smoothness of the dispersion function and tries to suppress the error along other directions (Izadian et al. 2019). By doing so, we not only managed to approximate the dispersion relation of a one-way wave equation, but also reduced the discretisation error of FD approximations. However, the FD equation (25) is just like the FD schemes of ABCs. Therefore, the proposed method is considered as a local semi-exact solution.

Numerical modelling
We compare the proposed boundary condition method with the second-order ABC (Clayton & Engquist 1977) and PML (Hustedt et al. 2004) (17). We use a Ricker wavelet with a dominant frequency of 20 Hz (Wang 2015) and maximum frequency of 60 Hz in the numerical examples. Note that a gain is applied to all the following figures to boost up small amplitudes of the reflections that could not be seen otherwise.
66 Figure 4. The energy of the propagating wavefield versus time of propagation. The proposed absorbing boundary condition performs better than the PML method with 10 nodes (PML-10).

Homogeneous media
The first example is a velocity model with 1.5 km × 1.5 km dimensions. A source wavelet is located at the centre of the model. To prevent numerical dispersion error be present in the modelling, we set the grid interval to be 5 m (six nodes per minimum wavelength). The thickness of the PML layer is set to be 10 nodes. The P-wave velocity is 2000 m s −1 . Figure 3 parts a-c show the real part of wavefields at 20 Hz, using the ABC, PML and the proposed boundary condition method, respectively. Three images were plotted with the same colour scale. The ABC result contains radial fringes associated to the reflections from boundaries (figure 3a), while the PML and the proposed method are free of these artificial reflections. Figure 3 parts d-f show the snapshots of the wavefields in time-domain at 0.5 s. The proposed method shows an acceptable performance for absorbing the boundary reflections.
To have a better quantitative evaluation of the performance of the methods, figure 4 compares the energy of wavefields (L 2 -norm of the wavefields) while propagating  in time, which again shows that the proposed method is performing better than the PML method.
The homogenous model test shows the efficiency of the proposed method in comparison with the ABC and PML methods. To evaluate the performance of the PML method with different thicknesses against the proposed absorbing condition, we performed another test similar to figures 3 and 4. We run the simulation with three PML layer thicknesses (10, 20 and 30 nodes), and show the results in figures 5 and 6. The results demonstrate that the proposed method performance is comparable with PML-30. In computationally expensive applications, padding boundaries of a model with extra 30 nodes requires more memory and extra calculations, which is not desirable. In contrast, our proposed method does not add new complexity to the problem.

A rectangular-shaped model
In the second example, we investigate performance of the proposed absorbing boundary method for wide angle reflections. For this purpose, we use a rectangular model with dimensions of 1.25 km × 3.25 km. Model velocity and simulation parameters are same as the previous example. The source is located at (1150 m, 1625 m), 100 m above the bottom edge of the model. Figure 7 parts a-c show the snapshots of wavefield at 0.6 s. The proposed method shows a better performance in suppressing unwanted reflections from the bottom boundary of the model.
A shot gather of the propagating wavefield at a depth of 1150 m captured reflections with incident angles of approximately 85 0 (figure 8). Following the primary wavefront, we see small reflections from the bottom edge. The ABC's result is the worst one as the reflection is obviously visible. PML's result (with 20 nodes) is showing the typical behaviour of this family, i.e. small low-frequency reflections (Gao et al. 2017). Finally, our proposed method indicates a higher absorbing performance as almost no tail is visible.

Anticline model
The third example is a layered model. The velocity increases with depth from 1500 to 2800 m s −1 (figure 9a). The grid spacing is 5 m. The source location is at the surface at x = 1375 m. In this simulation, we considered a free boundary condition for the top boundary and absorbing boundaries for other boundaries.
The snapshots of the simulated wavefield are shown for three different absorbing boundary methods in figure 9b-d. Artificial reflections are highlighted by black arrows in these figures and are more visible for the second-order ABC case. According to this figure, the proposed method could attenuate artificial reflections from all of boundaries effectively. Comparison of the results reveals that the proposed method has a better performance than the 10-point thickness PML in attenuating strong reflection from bottom boundary. Figure 10 shows shot gathers in which the receivers are located at top boundary for different absorbing boundary approaches. The proposed absorbing boundary method can attenuate boundary reflections effectively even better than a 10-point standard PML.
For better comparison, a trace located at x = 1150 m is extracted from the shot record ( figure 11). There is no significant difference between two results except the reflection from bottom of PML layer that is successfully suppressed by the proposed method.

Marmousi model
The last example deals with a highly complex velocity model, a part of the Marmousi model ( figure 12a). The grid interval is 8 m. By looking at the black arrows in the snapshots of the propagated wavefield in figure 12b-d, we can see a similar pattern of artificial boundary reflections from ABC and almost no reflection from PML20 and proposed method. In figure 13, the recorded data at the surface are shown. Black arrows in this figure highlight the boundary reflections.
6. Discussion on the computational efficiency As described, we used only one layer of neighbour grids in the FD schemes for both the computational domain's grids and boundary grids. Therefore, the structure of the impedance matrix will not be altered and, consequently, the computational cost will not be increased compared to the low-order ABCs. Figure 14 parts a-c show the structure of the impedance matrix of a 5 × 5 grid model for the conventional first, second-order and the proposed ABC methods at all edges. The blue points in figure 14c are the additional nonzero entries corresponding to the proposed in which all the blue points fall in the same three-band structure of the impedance matrix of the conventional ABCs. Hence, the computational cost of the proposed method is almost equal to a second-order ABC, which is almost the cheapest boundary condition available.  This very high performance of the proposed method will have a significant impact on the efficient frequency-domain FD modelling in sensitive imaging applications such as FWI where an accurate ABC is crucial. In a 3D application, a model with the size of 4 × 4 × 10 km, discretised with 20 m grid interval, will end up with 200 × 200 × 500 grids. Adding PML layers with thicknesses of 10, 20, and 30 grids will increase the number of grids by 20, 42 and 67%, respectively, whereas the proposed method, as a local absorbing boundary, increases the number of grids by a mere 1.9%. 71 Figure 14. Structure of the impedance matrix for a 5 × 5 grid model with a nine-point FD scheme for the interior domain grids and the first-order ABC (a), the second-order ABC (b) and the proposed method (c) as the boundary condition at all of edges.

Conclusions
We have proposed a local absorbing boundary condition for the seismic wave simulation using the FDFD solution of the acoustic wave equation. The observations from numerical analyses have confirmed the following: (1) The proposed method performs even better than the standard PML in a wide range of incident angles.
(2) The proposed method performs as well as the PML method having 30 nodes thickness with the difference that our boundary condition does not put the extra computational burden to the problem. (3) Despite being a local boundary condition, the performance of the suggested method was preserved even in a complex structured model such as the Marmousi model. (4) Furthermore, the local characteristics of this method makes it computationally efficient, which makes it an ideal choice for a broad application prospect of FWI.