Frequency-domain elastic wavefield simulation with hybrid absorbing boundary conditions

The frequency-domain seismic modeling has advantages over the time-domain modeling, including the efficient implementation of multiple sources and straightforward extension for adding attenuation factors. One of the most persistent challenges in the frequency domain as well as in the time domain is how to effectively suppress the unwanted seismic reflections from the truncated boundaries of the model. Here, we propose a 2D frequency-domain finite-difference wavefield simulation in elastic media with hybrid absorbing boundary conditions, which combine the perfectly matched layer (PML) boundary condition with the Clayton absorbing boundary conditions (first and second orders). The PML boundary condition is implemented in the damping zones of the model, while the Clayton absorbing boundary conditions are applied to the outer boundaries of the damping zones. To improve the absorbing performance of the hybrid absorbing boundary conditions in the frequency domain, we apply the complex coordinate stretching method to the spatial partial derivatives in the Clayton absorbing boundary conditions. To testify the validity of our proposed algorithm, we compare the calculated seismograms with an analytical solution. Numerical tests show the hybrid absorbing boundary condition (PML plus the stretched second-order Clayton absorbing condition) has the best absorbing performance over the other absorbing boundary conditions. In the model tests, we also successfully apply the complex coordinate stretching method to the free surface boundary condition when simulating seismic wave propagation in elastic media with a free surface.


Introduction
Frequency-domain seismic wavefield simulation has been developed over the last three decades (Marfurt 1984;Marfurt & Shin 1989;Pratt 1990;Pratt & Worthington 1990;Jo et al. 1996). Compared with seismic modeling in the time domain, frequency-domain seismic wavefield simulation has the advantages of efficient implementation of multiple sources and straightforward extension for adding attenuation factors. It also makes simulation of seismic wave propagation possible with only a few frequency components (Marfurt & Shin 1989). Numerous studies have been conducted in the frequency domain, such as seismic wavefield simulation (Song et al. 1995;Jo et al. 1996;Operto et al. 2007Operto et al. , 2009Zhao et al. 2017), full waveform inversion (Ben-Hadj-Ali et al. 2011;Tao & Sen 2013;Chi et al. 2014) and reverse time migration (Kim et al. 1996;Chung et al. 2012;. One of the most persistent problems in frequency-domain seismic modeling as well as time-domain modeling is how to effectively absorb artificial reflections that are caused by the truncated computational bounds of the model (Chen 2011). If the absorbing boundary condition does not work well in the frequency domain, there will exist serious artificial reflections in the wavefield solution of each frequency component. This further causes large errors in the time domain because the seismic record in the time domain can be regarded as the sum of all frequency components (Shin 1995).
Various absorbing boundary conditions have been applied to frequency-domain seismic modeling, including Clayton absorbing boundary conditions (Clayton & Engquist 1977) and the perfectly matched layer (PML) (Berenger 1994). However, each absorbing boundary condition has limitations and cannot be perfect (Tsynkov 1998).
Since Clayton & Engquist (1977) proposed the first-and second-order absorbing boundary conditions based on paraxial approximations of the scalar and elastic wave equations, they have been widely used in the seismic wavefield simulation in the time domain (Sheen 1993) and frequency domain (Pratt 1990;Song et al. 1995;Kim et al. 1996;Doyon & Giroux 2014;Moreira et al. 2014). The idea is very simple and only considers the outward-going wavefield after separating the inwardand outward-going wavefields. Clayton & Engquist (1977) also indicated that the Clayton absorbing boundary conditions are valid at both low and high frequencies. However, the Clayton absorbing boundary conditions only work well at normal incidence (Engquist & Majda 1977, 1979Cerjan et al. 1985;Renaut & Petersen 1989;Gao et al. 2017). In addition, for the special treatments at the four corner points, Clayton absorbing boundary conditions are derived by assuming a 45°incident angle, which may cause instability problems (Engquist & Majda 1979). Berenger (1994) introduced the PML absorbing boundary condition, which was theoretically proven to have a zero reflection coefficient for all angles of incidence without the discretization for electromagnetic wave equations in the time domain. To implement the PML boundary condition, we use a complex coordinate stretching method to the original coordinates by adding the linearly growing imaginary part, which makes the wave exponentially decay. After the PML boundary condition was applied to the seismic wavefield simulation in the time domain (Collino & Tsogka 2001), it was further used in frequency-domain seismic modeling (Operto et al. 2007;Liao et al. 2009;Li et al. 2015;Zhao et al. 2017). Although the PML boundary condition has the advantage of absorbing the evanescent field in discrete space, it suffers from spurious reflections at grazing incidence (Komatitsch & Martin 2007) and instability problems with the existence of surface waves (Festa et al. 2005). Additionally, the standard PML boundary condition encounters the instability problem on the interaction between the low-frequency evanescent waves and artificial layers (Festa et al. 2005;Gao et al. 2017).
To overcome the limitations of the existing boundary conditions, some hybrid boundary conditions have been developed and have shown great absorbing behaviors (Liu & Sen 2010;Ren & Liu 2013;Liu et al. 2017a). Ren & Liu (2013) extended the hybrid absorbing boundary condition proposed by Liu & Sen (2010) in the time domain to the frequency domain, and verified the ability of absorbing artificial reflections in different types of media. The main idea is to divide the computational domain into inner and transition areas. The performance of their proposed algorithm depends on how many transition areas are implemented in the model. The proposed algorithm requires high computational cost in the linear weighting process within the transition zones. Also, the weighting coefficient is hard to decide. Moreira et al. (2014) proposed a hybrid absorbing boundary condition that combines the PML and the first-order Clayton absorbing boundary condition for the acoustic wavefield simulation in the frequency domain. The hybrid absorbing boundary condition shows a good reduction in spurious reflections.
In this article, we implement hybrid absorbing boundary conditions to simulate frequency-domain elastic wave propagation to overcome the limitations of the Clayton (instability problem at the corner parts and poor performance for wave traveling at grazing incidence) and the PML (e.g. instability problems with the existence of surface wave and low-frequency waves) absorbing boundary conditions. We apply the traditional PML boundary condition in the damping areas of the model and the Clayton absorbing boundary conditions (first-and second-orders) at the outer boundaries of the damping areas. To improve the absorbing behaviors of the hybrid absorbing boundary conditions and the interaction between the Clayton and PML absorbing boundary conditions (not by simply 'adding' them together), we also apply the complex coordinate stretching method to the spatial partial derivative equations of the Clayton absorbing boundary conditions. Similarly, we apply the complex coordinate stretching to the free surface boundary condition in the upper left and right damping zones to improve the interaction between the free surface and the hybrid absorbing boundary conditions. Several numerical algorithms can be used to simulate seismic wave propagation in the frequency domain, such as the finite element method (Marfurt 1984;Zhao et al. 2017), pseudo-spectral method (Yan & Liu 2016), boundary element method (Chen & Zhou 1994), spectral element method (Komatitsch et al. 2000;Seriani & Oliveira 2007) and finite difference method (Moreira et al. 2014). Here, we choose the finite-difference method to carry out the frequency-domain seismic wavefield simulation because it provides a good balance between the numerical accuracy and computational cost. In addition, the Intel Pardiso Solver is applied to solve the large sparse matrix encountered in the frequency-domain seismic modeling (Liu et al. 2017b).

Frequency-domain elastic wave equations with the PML absorbing boundary condition
The main idea of the PML is to introduce complex coordinate stretching into a normal coordinate system, which adds the linearly growing imaginary part to exponentially decay the wavefield. For instance, the complex coordinate stretching along x direction can be expressed bỹx → 1 e x x , where e x is the coordinate stretching function andx is the stretched coordinate. By using the complex coordinate stretching method, the implementation of the PML absorbing boundary condition in the frequency-domain elastic wave equations can be expressed as (Liao et al. 2009): and We obtain and where x i and z i are the distances from the PML boundary in the horizontal and vertical directions, respectively; a 0 is the optimized parameter subject to the change of thickness, which is normally referred to 1.79 (Zeng et al. 2001), L PML is the thickness of the PML boundary and f 0 is the dominant frequency of the source. In this article, after our numerical tests, the best value of a 0 is 1.79 for the PML thickness of 30 grids.

First-and second-order Clayton absorbing boundary conditions
The first-and second-order absorbing boundary conditions proposed by Clayton & Engquist (1977) are, respectively, written as: where V p and V s are the P-and S-wave velocities, respectively, i is the complex imaginary part, is the angular frequency and ⃗ u is the displacement vector. The first-order Clayton absorbing boundary condition (equation (5)) has the advantage of time efficiency. However, it only works well for the wave propagation at normal incidence. Although the second-order Clayton absorbing boundary condition (equation (6)) performs better than the first-order Clayton absorbing boundary condition at other incidence angles, it may encounter some instability problem at corners of the model (Engquist & Majda 1979).
For the corner points of the model, the first-order Clayton absorbing boundary condition is rotated by 45º, which is expressed as: As we cannot rotate the second-order Clayton absorbing boundary conditions by 45º due to the limited number of points, we apply the same boundary condition indicated by equation (7) to the corner points of the model when carrying out seismic 692 wavefield simulation with the second-order Clayton absorbing boundary condition. In addition, the discretized forms of the first-and second-order Clayton absorbing boundary conditions also fit the finite-difference star stencil (Pratt 1990).
To improve the absorbing performance of the hybrid absorbing boundary conditions in our algorithm, we apply the complex coordinate stretching method to the first-and second-order Clayton absorbing boundary conditions, respectively. We obtain the following forms: and (i ) Similarly, we obtain the following form of the first-order Clayton absorbing boundary condition at the corner points of the model:

Free surface boundary condition
To model seismic wave propagation at the surface (air/solid interface), we apply the traction free surface boundary condition to the surface of the model (Vidale & Clayton 1986), which can be expressed as: and .
In the implementation of the free surface boundary condition, a fictitious layer is added above the actual free surface to allow the implementation of the finite differences of horizontal derivatives in the model (figure 1). As we do for the Clayton absorbing boundary conditions, the complex coordinate stretching method is also applied to the free surface boundary condition at the left and right boundaries of the free surface.

Numerical experiments
Here, two numerical tests are carried out to testify the performance of the proposed hybrid absorbing boundary conditions in the frequency-domain elastic wavefield simulation.

Case 1: homogeneous elastic model
The homogeneous elastic model (figure 2) is discretized by 300 × 300 grid points with grid space of 4 m. The P-wave velocity, S-wave velocity and density of the model are 3000 m s −1 , 2000 m s −1 and 2000 kg m −3 , respectively. The thickness of the PML boundary is 120 m (30 grid points). We use the same thicknesses for hybrid absorbing boundaries. The frequency range in the wavefield simulation is from 0.1 to 80.1 Hz with an increment of 1 Hz. A vertical force source (star) with a peak frequency of 20 Hz is located at the coordinates (0.60 km, −0.30 km). Two receivers (triangles) are located at coordinates (0.16 km, −0.16 km) and (1.00 km, −0.02 km), which are used to test the performance of absorbing boundary conditions at small and large incidence angles, respectively. In figure 2, the dashed lines represent the inner boundaries of the PML. Figure 3 shows the 30.1 Hz (left panel) and 70.1 Hz (right panel) real-part solutions of the vertical component using different boundary conditions. One can observe the strong artificial reflections in the wavefields using the first-order (figure 3a and b) and second-order (figure 3c and d) Clayton absorbing boundary conditions. However, we do not observe the distinct artificial reflections in the wavefields by using the PML boundary condition (figure 3e and f), the hybrid absorbing boundary condition (PML plus the first-order Clayton boundary condition) (figure 3g and h) and the hybrid absorbing boundary condition (PML plus the second-order Clayton boundary condition) ( figure 3i and j).
To further test the absorbing performance of our proposed hybrid absorbing boundary conditions, we compare seismograms of the vertical component for receivers R1 (figure 4) and R2 (figure 5), respectively, with an analytical solution (Aki & Richards 2002). We obtain the seismic wavefields in the time domain that are converted from the frequency domain by using an inverse Fourier transform. We can observe that the seismograms (figures 4f and 5f) calculated by the hybrid absorbing boundary condition (PML plus the stretched second-order Clayton absorbing boundary condition) have the best agreement with the analytical solution. The cross-correlation parameters (XCORR) of the calculated seismograms with the analytical solution is used to verify the performance of the proposed algorithm quantitatively. We can see that the cross-correlation values in regard to the hybrid boundary condition of the PML and stretched second-order Clayton absorbing boundary condition are the highest. Figure 6 represents the wavefield energy decay curves (up to 1 s) by using different boundary conditions. After the Pwave reaches the absorbing boundaries of the model (theoretically at 0.06 s), the energy carried by seismic wave is gradually absorbed. After approximately 0.46 s, theoretically the energy remained in the model should be caused by the artificial reflections because both P-and S-waves have left the main computational domain. One can observe that the energy calculated with the hybrid absorbing boundary condition (PML plus the stretched second-order Clayton absorbing boundary condition) is absorbed faster than the energy calculated with other absorbing boundary conditions.

Case 2: homogeneous elastic model with a free surface
To test the performance of the hybrid absorbing boundary conditions with the existence of a surface wave, we built a homogeneous elastic model with a free surface ( figure 7). The model has a size of 0.9 × 0.45 km, which is discretized by 300 × 150 grid points with a grid space of 3 m. The frequency range in the wavefield simulation is from 0.1 to 79.85 Hz with an increment 694  (Vidale & Clayton 1986) is only applied at the surface. The star denotes the horizontal force source with a peak frequency of 20 Hz located at the coordinates (0.15 km, −0.006 km), and the triangle denotes a receiver located at the coordinates (0.16 km, −0.08 km). one can easily observe that there is no spurious reflection from the truncated computational domain. We can see a significant reduction in the artificial seismic reflections in figures 8e and 8f. In addition, figure 8 shows that our proposed hybrid absorbing boundary condition can handle the free surface at the top and have good absorbing effects.
To better observe the performance of the proposed hybrid absorbing boundary condition, we also show seismic profiles (figure 9) calculated with different boundary conditions in this study. Figure 9 shows seismic profiles calculated with a second-order Clayton absorbing boundary condition (figure 9a), PML absorbing boundary condition (figure 9b) and hybrid boundary condition (PML plus the stretched second-order Clayton absorbing boundary condition). One can also observe that our algorithm can handle the free surface, and the proposed hybrid absorbing boundary condition has better absorption performance than the second-order Clayton absorbing boundary condition and the PML absorbing boundary condition.  We also calculate the seismograms of the horizontal and vertical components for the receiver at (0.75 km, −0.03 km). For simplicity, figure 10 only shows the seismograms calculated with the second-order Clayton absorbing boundary condition ( figure 10a and b), the PML (figure 10c and d) and the corresponding hybrid absorbing boundary condition (figure 10e and f), respectively. We can find strong spurious reflections in the seismograms (figure 10a and b) calculated with the second-order Clayton absorbing boundary condition. Although the PML boundary condition (figure 10c and d) works better, the hybrid absorbing boundary condition (figure 10e and f) provides the best seismogram compared with other absorbing boundary conditions. Figure 11 represents the wavefield energy decay curves (up to 3 s) by using different boundary conditions. After the P-wave reaches the absorbing boundaries of the model, energy carried by the seismic wave is gradually absorbed. After the Rayleigh wave reaches the truncated computational boundaries, the energy remaining in the domain is caused by the artificial reflections from the boundaries. In addition, one can observe that the energy calculated with the hybrid absorbing boundary condition (PML plus the stretched second-order Clayton absorbing boundary condition) is absorbed faster than the energy calculated with other absorbing boundary conditions.

Discussion
In this section, we first discuss the performance of the PML absorbing boundary condition with different thicknesses and different damping coefficients. Then, we discuss the efficiency of the proposed algorithm in regard to grazing incidence and complex media: modified Marmousi-ii model (Martin et al. 2002). The computational costs of the computing codes with various boundary conditions are discussed at the end of this part.   Figure 12 shows the wavefield energy decay curve for different damping coefficients based on the homogeneous elastic model (case 1). We can observe that the damping coefficient of 0 = 1.79 can provide the best energy absorption.

Influence of thickness and damping coefficient on the performance of the PML boundary condition
We further discuss the performance of the PML boundary condition with different thicknesses under the condition of the optimal damping coefficient. Figure 13 shows the horizontal component of seismic records along the surface calculated with 10-grid (a), 20-grid (b) and 30-grid (c) thicknesses of the PML boundary condition and 20-grid (d) thickness of the hybrid boundary condition (PML plus the stretched second-order Clayton absorbing boundary condition). In figure 13, we can see 699 Figure 10. Seismograms of the horizontal (left) and vertical (right) components recorded by the receiver at (0.75 km, −0.03 km) calculated using the second-order Clayton boundary condition (a and b), the PML absorbing boundary condition (c and d) and the hybrid absorbing boundary condition (PML plus the stretched second-order Clayton boundary condition) (e and f).  in figure 14. One can see that the absorbing performance of the hybrid boundary condition (PML plus the stretched secondorder Clayton absorbing boundary condition) with 20-grid thickness provides an excellent result compared to the PML boundary condition with 30-grid thickness.
For better absorbing performance of the hybrid absorbing boundary condition, we use the 20-grid thickness and damping coefficient of the PML boundary condition in the following discussions.

Grazing incidence
Here, we also choose the same homogeneous elastic model (case 1). The only difference is that the source of Ricker wavelet with the dominant frequency of 20 Hz and time delay of 0.2 s is located at (1.00 km, −0.20 km). For simplicity, we only show seismograms (figure 15) calculated with the second-order Clayton absorbing boundary condition, the PML boundary condition and the hybrid boundary condition (PML plus the stretched second-order Clayton absorbing boundary condition), which are recorded by the receiver located at (0.20 km, −0.24 km). In figure 15, we can see that the second-order Clayton absorbing boundary condition does not work well to suppress the artificial reflections in the case of the grazing incidence. The hybrid boundary condition can perform good absorption of the unwanted reflections at grazing incidence.

Complex elastic model: modified Marmousi-ii
We also test the performance of the proposed hybrid absorbing boundary condition in the complex elastic model: modified Marmousi-ii (Martin et al. 2002). The model has a size of 0.5 × 0.5 km, which is composed of 500 × 500 grids with a grid size of 1 m. Figure 16 illustrates the images of P-and S-wave velocities and density of the modified Marmousi-ii model.  second-order Clayton absorbing boundary condition) provides the best absorption of the artificial reflections in the computational domain (figure 17h) compared with other boundary conditions. Figure 18 shows the differences between figures 17a and 17c (figure 18a), and between figures 17e and 17g (figure 18b), which indicate the great differences in absorption performance among these four absorbing boundary conditions.
We also show the wavefield energy decay curves using different absorbing boundary conditions in figure 19. One can observe that the second-order Clayton absorbing boundary condition (dark dashed line) and the PML boundary condition (gray solid line) have a similar absorbing performance. Although both hybrid boundary conditions work well, the hybrid boundary condition (PML plus the stretched second-order Clayton absorbing boundary condition, dark solid line) has a better absorbing performance than the hybrid boundary condition (PML plus the unstretched second-order Clayton absorbing boundary condition, gray dashed line).    Table 1 shows the consuming times of computing codes with different boundary conditions in two cases used in the part of numerical experiments. All computing codes are running by the same Dell Workstation (Dell Precision T7500, 64-bit Dual Processor). From the table, we can see that the proposed algorithms of hybrid absorbing boundary conditions run with shorter consuming times than the algorithms of the classical PML boundary condition, and the hybrid absorbing boundary condition (PML plus the stretched second-order Clayton absorbing boundary condition) provides a very good balance between the computational cost and absorbing efficiency.

Conclusions
In this article, the performances of the hybrid absorbing boundary conditions, which respectively combine the first-and second-order Clayton absorbing boundary conditions with the PML boundary condition, were successfully investigated in a frequency-domain elastic wavefield simulation. Unlike previous studies, here we adopted the complex coordinate stretching method to the Clayton absorbing boundary conditions and free surface boundary condition. Numerical experiments show that the hybrid absorbing boundary condition (PML plus the second-order Clayton absorbing boundary condition) has the best performance in reducing spurious reflections when compared with other absorbing boundary conditions.