Theoretical and quantitative evaluation of hybrid PMLABCs for seismic wave simulation

Yuhang Wang and Wei Zhang,?

1 Shenzhen Key Laboratory of Deep Offshore Oil and Gas Exploration Technology, Southern University of Science and Technology, Shenzhen 518055, China

2 Southern Marinea Science and Engineering Guangdong Laboratory (Guangzhou), Guangzhou 511458, China

3 Department of Earth and Space Sciences, Southern University of Science and Technology, Shenzhen 518055, China

ABSTRACT A good artificial boundary treatment in a seismic wave grid-based numerical simulation can reduce the size of the computational region and increase the computational efficiency, which is becoming increasingly important for seismic migration and waveform inversion tasks requiring hundreds or thousands of simulations. Two artificial boundary techniques are commonly used:perfectly matched layers (PMLs), which exhibit the excellent absorption performance but impose a greater computational burden by using finite layers to gradually reduce wave amplitudes; and absorbing boundary conditions (ABCs), which have the high computational efficiency but are less effective in absorption because they employ the one-way wave equation at the exterior boundary. Naturally, PMLs have been combined with ABCs to reduce the number of PMLs, thus improving the computational efficiency; many studies have proposed such hybrid PMLs. Depending on the equations from which the ABCs are derived, there are two hybrid PML variants: the PML+unstretched ABC (UABC), in which the ABC is derived from a physical equation; or the PML+stretched ABC (SABC), in which the ABC is derived from the PML equation. Even though all the previous studies concluded that hybrid PMLs can improve the absorption performance, none of them quantified how many PMLs can be removed by combining the PML with the ABC compared with the pure PML. In this paper, we systematically study the absorption performance of the two hybrid PML variants. We develop a method to distinguish the artificial reflections from the PML-interior interface and those caused by the PML exterior boundary to accurately approximate the additional absorption achieved by using the UABC and the SABC. The reflection coefficients based on a theoretical derivation and numerical tests both show that the UABC amplifies most reflections and is not recommended in any situation; conversely, the SABC can always diminish reflections, but the additional absorption achieved by the SABC is relatively poor and cannot effectively reduce the number of PMLs. In contrast, we find that simply increasing the damping parameter improves absorption better than the PML+SABC. Our results show that the improvement in absorption achieved by combining the PML with either the SABC or the UABC is not better than that obtained by simply adjusting the damping profile of the PML; thus, combining the PML with the ABC is not recommended in practice.

Keywords: absorbing boundary condition; perfectly matched layer; Higdon boundary; hybrid; finite difference.

1. Introduction

Seismic migration and waveform inversion require considerable amounts of computation time and memory for artificial boundary treatments to eliminate spurious reflections, especially within large-scale, complex 3D models (Gao YJ et al., 2017). Two artificial boundary techniques are commonly used: absorbing boundary layers,mainly referring to perfectly matched layers (PMLs),which can gradually absorb waves in finite layers so that the ultimate reflection amplitude is acceptable, and absorbing boundary conditions (ABCs), which allow waves to propagate only outward at the boundary by using the one-way wave equation.

Bérenger (1994) proposed the standard PML for Maxwell’s equations that can theoretically eliminate the reflections of incident waves at all angles and all frequencies. Subsequently, the concept of the PML was quickly extended to seismic wave modeling, and considerable progress was soon achieved (Chew and Liu QH, 1996; Hastings et al., 1996; Collino and Tsogka,2001). However, strong spurious reflections continued to appear for near-grazing incident waves, low-frequency waves and evanescent waves (Festa et al., 2005;Komatitsch and Martin, 2007; Drossaert and Giannopoulos, 2007b). To mitigate these reflections,Kuzuoglu and Mittra (1996) proposed the complex frequency-shifted PML (CFS-PML), which was originally implemented in a split-field form (Gedney, 1998), after which some unsplit implementations were proposed by using convolutional terms (C-PML, Roden and Gedney,2000; Drossaert and Giannopoulos, 2007a; Komatitsch and Martin, 2007), integral terms (RI-PML, Zeng YQ and Liu QH, 2004; Drossaert and Giannopoulos, 2007b), the matched Z-transform (MZT CFS-PML, Shi RQ et al.,2012), and auxiliary differential equations (ADE CFSPML, Wang LN and Liang CH, 2006; Martin et al., 2010;Zhang W and Shen Y, 2010). Furthermore, to overcome the instabilities arising in some anisotropic media and to absorb waves over broad frequency bands, the multiaxial PML (MPML, Meza-Fajardo and Papageorgiou, 2008;Zhang ZG et al., 2014; Zhao ZC et al., 2019) and the highorder CFS-PML (Correia and Jin JM, 2005; Martin et al.,2010; Feng HK et al., 2017) were also adopted.

Generally, the PML can efficiently absorb waves with broad incident angles; however, this increases the size of the computational domain and thus demands considerable amounts of computation time and resources. To improve the computational efficiency, some authors have adjusted the damping profile of the PML to improve the absorption efficiency so that fewer PMLs can be used. The optimal damping profile shapes and the optimal parameters of PMLs have also been widely studied (Collino and Tsogka,2001; Appel? and Kreiss, 2006; Modave et al., 2010;Kucukcoban and Kallivokas, 2011; Yang XN and Zhang W, 2020). Another feasible scheme is to truncate the PML by ABCs, which are easy to understand and implement,impose a low computational demand, and achieve good absorption.

After many years of development, many kinds of classic ABCs are available, such as the Clayton-Engquist boundary (Clayton and Engquist, 1977), multi-transmitting formula (Liao ZF et al., 1984), and Higdon boundary(Higdon, 1986, 1987, 1992, 1994); other examples including the potential decomposition method (Randall,1988) and the arbitrarily wide-angle wave equations(Guddati and Heidari, 2005; Heidari and Guddati, 2006).In addition, high-order ABCs based on auxiliary variables were developed to easily absorb waves with a broad range of incident angles (Hagstrom and Warburton, 2004;Hagstrom et al., 2008). In recent years, some novel implementations of ABCs have been presented, such as multiple absorbing surfaces (MASs, Sudiarta, 2003) and the re-radiating boundary condition (rRBC, Diaz and Scherbatko, 2004, 2005), which were generalized by Bérenger (2007) as the Huygens absorbing boundary condition (HABC). Among all the existing ABCs, the Higdon boundary is one of the most widely used ABCs in practice because it is easy to extend to high orders and does not require special treatments for corners.

Some scholars have sought to combine PMLs with ABCs, and good results have been obtained. These combinations can be divided into two kinds based on the equations from which the ABC is derived: unstretched ABCs (UABCs), derived from the physical equation(Pekel and Mittra, 1995; Kantartzis and Tsiboukis, 1997;Xiao Y and Lu YL, 2001; Moreira et al., 2014; Darvish, et al., 2018) and stretched ABCs (SABCs), derived from the PML equation (Jin JM and Chew, 1996; Petropoulos,1998; Hu JL et al., 2018; Zhao ZC et al., 2019). Moreira et al. (2014) combined the PML with the unstretched Clayton ABC in acoustic wave modeling and reported that the hybrid scheme achieved better performance. Darvish et al.(2018) combined the PML with the unstretched Higdon ABC in electromagnetic wave modeling and determined that the hybrid method could lead to an approximate enhancement in absorption of 10 dB. Hu JL et al. (2018)combined the PML with the stretched Higdon ABC in acoustic modeling and also indicated that the hybrid boundary could absorb waves more efficiently. Zhao ZC et al. (2019) combined the PML with the stretched Clayton ABC in elastic wave modeling and concluded that the hybrid boundary provided a good balance between the computational cost and the absorption effectiveness.

Nevertheless, even though all the previous studies concluded that hybrid PMLs can improve the absorption performance, none of them quantified how many PMLs can be removed by combining the PML with the ABC compared with the pure PML. In this paper, we systematically compare the absorption improvement of the hybrid PMLs in combination with the first-order unstretched Higdon ABC (PML+UABC) and the firstorder stretched Higdon ABC (PML+SABC). We derive the theoretical reflection coefficients in continuous space and semidiscrete space to investigate whether the PML+UABC and PML+SABC can work from a theoretical perspective. We also develop a method to separate the artificial reflections caused by the PMLinterior interface and those caused by the PML exterior boundary to accurately reflect the additional absorption achieved by using the UABC and the SABC from the numerical results. Moreover, we compare the results with the absorption improvement obtained by simply adjusting the damping profile of the PML. Our results ultimately show that the absorption improvement achieved by combining the PML with either the SABC or the UABC is not better than that obtained by simply adjusting the damping profile of the PML; thus, the combination of the PML with the ABC is not recommended in practice.

2. Derivations of the PML, UABC and SABC

In this section, we briefly describe the implementation of the ADE CFS-PML and then derive the equations of the UABC and the SABC for the Higdon ABC.

2.1. ADE CFS-PML

We consider the following first-order velocity-stress equation:

wherevis the particle velocity vector, τ is the secondorder stress tensor, the superscript T is the transpose operator, ρ is the density of the medium, andCis the fourth-order stiffness tensor.

Based on the concept of complex coordinate stretching(Chew and Weedon, 1994), the equations in the PML domain can be easily derived by simply replacing the coordinatexwith the complex stretched coordinatex?. For simplicity, we define the start pointx=0 for the PML and discuss only the +xdirection here:

wheresxis the complex stretching function, ξ is the integral variable, αxis the frequency-shifted function, βxis the scaling function, anddxis the damping function that exponentially reduces the wave amplitude. Specifically,when αx=0 and βx=1, the CFS-PML degenerates into the standard PML.

In this paper, we implement the CFS-PML by using ADEs, which simplify the discretization process by using the same numerical scheme in the original equations and auxiliary equations. We take the expression of thevxcomponent of Equation (1) in thexdirection as an example:

whereωis the circular frequency, and ?· represents the temporal Fourier transform of the variable. We decouplesxinto two parts:

Referring to Zhang W and Shen Y (2010), Equation (3)can be written as:

2.2. Unstretched Higdon ABC

We briefly derive the unstretched Higdon ABC in continuous space, though the Higdon operator was originally defined in discrete space (Higdon, 1986). We consider the two-dimensional scalar acoustic wave equation:

whereuis the wavefield.

Applying temporal and spatial Fourier transforms to Equation (6) yields the following dispersion relation:

wherecis the wave velocity, andkxandkzare the spatial wavenumbers.

Thus, we have:

Assuming thatckz/ω is small, we expand the square root operator via a Taylor expansion and take the firstorder terms about the - sign:

Next, we apply an inverse Fourier transform to Equation (9):

To achieve better absorption at a certain angle α,Equation (10) can be extended as:

Finally, we can directly obtain the high-order unstretched Higdon ABC by multiplying several Higdon operators:

where αiandciare exactly compatible with the outward flux of energy via plane waves at any incident angle αiand any wave velocityci.

Although the derivation of the unstretched Higdon ABC above originates from the acoustic wave equations,this UABC can effectively handle elastic waves (Higdon,1991).

2.3. Stretched Higdon ABC

We follow the same procedures as in the last subsection to derive the stretched Higdon ABC by using the dispersion relation, which is different from the derivation by the PML plane wave solution (Hu JL et al.,2018) and the direct extension by applying the stretching function to the unstretched ABC (Zhao ZC et al., 2019).We consider the following two-dimensional scalar acoustic wave equation in the PML:

Similarly, we apply the spatial Fourier transform to Equation (13):

whereu=Fxz[u?] and Fxzrepresents the spatial Fourier transform.

To obtain the dispersion relation, we address the first term on the right-hand side of Equation (14), which can be simplified as:

where * represents convolution.

We assume thatsxis a constant aboutxand is denoteds:

Then, Equation (15) can be simplified as:

Thus, we have the dispersion relation for Equation(13):

The above equation can be written as:

Assuming thatckz/ω is small, we expand the square root operator via a Taylor expansion and take the firstorder terms about the - sign:

Similarly, the generalized stretched Higdon ABC can be expressed as:

where αiandciare exactly compatible with the outward flux of energy via plane waves at any incident angle αiand any wave velocityciin the PML.

Referring to the solution process of the ADE CFSPML, we can easily obtain the updating formulas of the first-order stretched Higdon ABC traveling in the +xdirection:

For α,β,din Equation (22), we adopt the value at the PML exterior boundary. Forcin Equation (22), we choose the compressional-wave velocity for simplicity, which can also help absorb shear waves (Higdon, 1991).

For simplicity, only the first-order Higdon ABC is discussed in this paper. For convenience, hereafter, the UABC represents the first-order unstretched Higdon ABC,and the SABC represents the first-order stretched Higdon ABC.

3. Origin of the reflections and derivations of the reflection coefficients

Before studying the effects of the combinations of the PML with the UABC or the SABC, it is necessary to introduce the origin of spurious reflections, as this information provides useful guidance in the next section.

Nissen and Kreiss (2011) divided the errors of the discretized PML problem in a truncated domain into three types:

1. The discretization error ε0, which comes from the equation in the interior domain.

2. The modeling error ε1, also called the PML exterior truncation error, which comes from solving the continuous PML equation with a finite PML width.

3. Numerical reflections ε2, which come from the discretized PML equation. Specifically, these reflections can be divided into two parts: reflections from the PMLinterior interface and reflections from the increasing absorption function.

In practice, ε1and ε2are the main reflections in most cases. A small damping parameter leads to large modeling errors due to insufficient absorption (Collino and Monk,1998; Collino and Tsogka, 2001). Conversely, a large damping parameter leads to large spurious reflections from the PML-interior interface due to the inadequate representation of the PML (Bermúdez et al., 2007; Skelton et al., 2007; Nissen and Kreiss, 2011) and to reflections inside the PML coming from the rapidly increasing damping profile (Kucukcoban and Kallivokas, 2011;Nissen and Kreiss, 2011; Yang XN and Zhang W, 2020).

In this paper, we assume that ε0is negligible and focus only on ε1and ε2. In numerical tests, the artificial reflections due to ε1and ε2are usually very close and sometimes even overlap each other. Because the effect of either the UABC or the SABC is to be reflected by the value of ε1, we separate ε1and ε2from the total reflection error.

3.1. Reflection coefficients for the continuous model

Here, we derive the reflection coefficients in continuous space, which can provide an intuitive understanding of whether the UABC and the SABC can work.

For simplicity, we assume that waves propagate in the(x,z) plane and reflect from the right-hand boundary(x=L). We consider the following plane wave solution including an incident wave and a reflected wave in the frequency domain:

wherecis the velocity of the traveling wave, θ is the azimuth angle,Ris the reflection coefficient, andrepresents the complex coordinate, defined by:

Then, we can derive the reflection coefficients for different boundaries.

For the first-order Higdon ABC without a PML, we set=x,x≥0. Equation (23) should satisfy boundary condition Equation (12) withn=0 in the frequency domain atx=L:

We substitute Equation (23) into Equation (25):

For the CFS-PML+Dirichlet boundary, Equation (23)should satisfy=0 atx=L:

For the CFS-PML+UABC, Equation (23) should also satisfy boundary condition Equation (25). Then, we substitute Equation (23) into Equation (25):

where

andRUABCrepresents the reflection coefficients caused by the UABC only.

Now, we consider a simple example where cosα=cosθ=0 and αL=0,βL=1:

Notably, |RUABC|is not zero in this case. Moreover,dLincreases until it reaches one. Thus, the UABC tends to fail in practice because the reflection amplitude is theoretically nonzero (even close to one).

For the CFS-PML+SABC, Equation (23) should satisfy boundary condition Equation (21) withn=0 atx=L:

We substitute Equation (23) into Equation (31):

Notably,RSABCequalsRABC; thus, the PML+SABC’s reflection coefficient is the product of the PML’s reflection coefficient and the ABC’s reflection coefficient. Now, we consider the same example with cosα=cosθ=0 and αL=0,βL=1; here, |RABC| is clearly zero. Thus, in theory,the SABC can work well outside of the PML domain,similar to the ABC outside of the physical (inner) domain.

3.2. Reflection coefficients for the semidiscrete model

We next derive the reflection coefficients in semidiscrete space to help us quantitatively evaluate the absorption performance of the UABC and SABC.

Let us consider the spatial discretization of the plane wave solution in the (x,z), plane for the PML model with a finite width. The semidiscrete fieldtakes its value at the discrete pointswhereiandjare the spatial indices and Δxand Δzare the spatial steps.

A constant absorption functionsi=sis used in the PML domain (0 ≤i≤N).

The semidiscrete plane wave solution ofcan be written as:

wherekixandkzprepresent the discrete wavenumbers of the+xdirection in the inner domain and the PML domain,respectively, and are equal to ω cosθ/cwhen the dispersion is negligible in the numerical tests.Rrepresents the total reflection coefficients of the PML and its exterior boundary.

For the CFS-PML+SABC, we set c osα=0 in Equation(31) and use the first-order backward difference operator to calculate the spatial derivative:

Substituting Equation (33) into Equation (34), we can obtain:

For the CFS-PML+UABC, we lets=1 in Equation(34), yielding:

For the CFS-PML+Dirichlet boundary, we let=0 in Equation (33), yielding:

BecauseRin Equation (33) represents the total reflection coefficients of the PML and its exterior boundary, to obtain the PML exterior boundary’s reflection coefficient, we define:

whereRUABCandRSABCrepresent the reflection coefficients of the UABC and SABC, respectively.

Fi nally, the moduli ofRUABCandRSABCcan be written as:

where

and θ is the azimuth angle of the plane wave.

Equations (40) and (41) are used to verify our numerical results in the next section.

4. Numerical experiments

4.1. Parameter options for the PML

Before conducting the numerical tests, we briefly summarize the options for the PML parameters. In this paper, for simplicity, we discuss the absorption performance of only the standard PML+UABC and the standard PML+SABC.

Fordx, we adopt the commonly usedp-order polynomial profile in the PML (Gedney, 1998; Roden and Gedney, 2000):

whered0is the maximum value ofdxat the exterior boundary of the PML,Lis the width of the PML, andxis the distance from the computation point to the PMLinterior interface.

Skelton et al. (2007) stated thatpdshould be ≥ 2,which can help avoid reflections from the PML-interior interface without applying any special treatment to the derivatives at the interface. Here, we use the widely used order numberpd=2.

Collino and Tsogka (2001) presentedd0through a theoretical reflection coefficient relation based on normally incident waves:

wherecprepresents the compressional-wave velocity;Rcan be given asR5=0.01,R10=0.001 andR20=0.0001; and the superscript ofRdenotes the number of PMLs.

Here, we adoptRas defined by Zhang W and Shen Y(2010) and use the variableR0to replace the constant:

whereR0=3 represents the option of Zhang W and Shen Y (2010), who approximated the relationship betweenRand the PML width developed by Collino and Tsogka(2001) as a logarithmic equation andR0=4 represents the option of Yang XN and Zhang W (2020), who modified the option of Zhang W and Shen Y (2010) based on many numerical tests.

In this paper, we mainly discuss two different options ofR0(∝d0) by choosingR0=3 orR0=4.

4.2. Origin of the reflections

Nissen and Kreiss (2011) used a constant damping profile to study the reflections from the PML-interior interface and ignored the reflections from the increasing damping profile. Here, we adopt the quadratic damping profile in our tests, which is commonly used in practice.Next, we determine the origins of the reflections by performing numerical tests, which can help us avoid possible interference when we study the absorption performance of the UABC and SABC.

The medium has a compressional-wave velocityvP=3000m/s and a density ρ=2000kg/m3. The grid spacing is 50 m, the total number of grids is 801, and the size of the model is 40 km. The source is added to the velocity and applied at 20 km. Thirty PMLs with the Dirichlet boundary are used to distinguish different reflections. The standard PML (α =0,β=1) andR0=3(d0=32.2,R=2.6×10-5) are used by default. The quadratic damping profile depicted by Equation (43) is used if not otherwise specified.

The source time function is a Ricker wavelet with a center frequency of 2 Hz and a delay of 0.5 s. Gaussian spatial smoothing is used for the source with a time interval of 0.01 s and a total of 1 000 time steps. Here, we use a fourth-order staggered finite difference (FD) operator to calculate the first spatial derivatives inside the inner(physical) domain and gradually lower the order at the boundary to the second order for the ease of implementing ABCs. The fourth-order Runge-Kutta scheme is adopted for the time-marching method.

Figure 1 shows a snapshot ofvxfor incident waves.Figure 2 shows a snapshot ofvxfor the reflections.Reflections 1 and 2 are clearly seen; next, we determine their origin by performing numerical tests. Considering the symmetry of the model, we take the right part (x≥400) for further analysis.

Figure 1. Snapshot of vx for incident waves. R0 = 3, the dashed line represents the PML-interior interface, and the horizontal arrows represent the traveling direction of the incident waves.

Figure 2. Snapshot of vx for reflections. R0 = 3, the dashed line represents the PML-interior interface, and the horizontal arrows represent the traveling direction of the reflections. We call the first reflection at x = 620 reflection 1 and the second reflection at x = 680 reflection 2.

Figure 3. Damping profile of the whole PML. R0 = 3, and the dashed line represents the PML-interior interface (x = 771).

Figure 4. Magnified view of the damping profile near the PML-interior interface. R0 = 3, and the dashed line represents the interface (x = 771).

First, we briefly verify that reflection 1 is caused by the PML-interior interface. Figure 3 shows the damping profile of the whole PML. Figure 4 shows a magnified view of the damping profile near the PML-interior interface (x= 771). We adopt the stain algorithm (Chen B and Jia XF, 2014; Chen B et al., 2016; Li QH and Jia XF,2017) to investigate where the artificial reflections originate. The stain algorithm changes the model’s parameters at target points (stained structures) to be complex and then simulates a stained wavefield activated by these stained structures as the secondary sources when the normal wavefield passes them. The stained wavefield can effectively indicate the kinematic information of waves that are scattered or reflected from the stained structures. In this test, we stain 5 points one by one from 771 to 773 at an interval of one-half of a grid point (for the different wavefield components at the staggered grid).Figure 5 shows a comparison between reflection 1 and the stained wavefields originating from different stained points. Because the amplitudes of the stained wavefields is not amplitude preserving, the waveforms are normalized in Figure 1. We also take the spatial derivative of the stained wavefield to fit the wavelet of reflection 1. It can be seen that only the stained wavefield fromx= 771 (the location of the PML-interior interface) matches reflection 1, which indicates that reflection 1 is caused by the PML-interior interface.

Figure 5. Comparison between reflection 1 and the stained wavefields originating from different stained points. R0 = 3, the solid line represents the original unstained wavefield, and the dashed lines represent the processed stained wavefields originating from x = 770, 770.5, 771, 771.5, and 772 from left to right.

Second, we verify that reflection 2 is caused by the PML exterior boundary through a simple test. Figure 6 shows a comparison of the normalized reflection caused only by the Dirichlet boundary with that caused by the PML+Dirichlet boundary. Reflection 2 completely overlaps with the reflection caused only by the Dirichlet boundary, which means that reflection 2 is caused by the PML exterior boundary.

Figure 6. Comparison of the vx for the normalized reflections caused only by the Dirichlet boundary and that caused by the PML+Dirichlet boundary with R0 = 3.

These two tests verify two different reflections:reflection 1 is caused by the PML-interior interface,whereas reflection 2 is caused by the PML exterior boundary. Notably, no reflections are visible from the increasing absorption function.

Next, to test whether any reflections originate from the rapidly increasing damping profile, we choose a relatively largeR0=7 (d0=60.3, representing the sharpest damping profile among our tests) and repeat the procedures above.Figure 7 shows a snapshot ofvxfor the reflection withR0=7. Owing to the considerable absorption achieved by the rapidly increasing damping profile of the PML,reflection 2 is invisible, whereas reflection 1 is dominant.Figure 8 shows a comparison between reflection 1 and the stained wavefields originating from different stained points, and Figure 9 shows a comparison of the normalized reflection caused only by the Dirichlet boundary with that caused by the PML+Dirichlet boundary withR0=7. The conclusion is the same as that withR0=3. Thus, there is no need to consider the effects of reflections from the increasing damping profile in our tests. Notably, for largeR0, reflection 1 dominates, and thus, its effects must be removed when we study the absorption performance of the UABC and the SABC, which are related to reflection 2.

Figure 7. Snapshot of vx for reflections. R0 = 7, the dashed line represents the PML-interior interface, and the horizontal arrows represent the traveling direction of the reflections.

Figure 8. Comparison between reflection 1 and the stained wavefields originating from different stained points. R0 = 7, the solid line represents the original unstained wavefield, and the dashed lines represent the processed stained wavefields originating from x = 770, 770.5, 771, 771.5, and 772 from left to right.

Figure 9. Comparison of the vx for the normalized reflections caused only by the Dirichlet boundary and that caused by the PML+Dirichlet boundary with R0 = 7.

Figure 10. Comparison among the reflections of vx for the PML+Dirichlet boundary, PML+UABC and PML+SABC with R0 = 3; 30 PMLs are used. The reflections at approximately x =630 and x = 690 are due to ε2 and ε1, respectively.

4.3. Absorption performance of the PML+UABC and PML+SABC in 1D tests

In this subsection, we compare the absorption performance of the PML+UABC and PML+SABC by performing 1D numerical tests. The parameters are the same as in the last subsection.

For convenience, we define the reflection coefficients of the UABC and the SABC in numerical tests as:

Figure 10 shows a comparison among the reflections ofvxfor the PML+Dirichlet boundary, PML+UABC and PML+SABC. Reflections 1 and 2 are totally separate. The maximum amplitudes of reflection 2 for the PML+Dirichlet boundary, PML+UABC and PML+SABC are 113, 140,and 22, respectively. The PML+UABC amplifies the reflection, withreaching 1.24, while the PML+SABC weakens the reflection, andis 0.19.Notably, Equation (32) indicates thatRSABCequalsRABCin continuous space; however, in discrete space,differs too much from the reflection amplitude of the ABC(without the PML), which is generally less than 0.1.Thereafter, we derive the semidiscrete reflection coefficients of the UABC and the SABC [Equations (40) and(41), respectively] and find that the PML parameters have important impacts on these coefficients in discrete space.

Figure 11 compares the |R| derived from the semidiscrete plane wave solution with the |R| calculated by the numerical tests for the UABC and the SABC as a function ofd0. The center frequency 2 Hz is used to calculate the semidiscrete |R| in Equations (40) and (41). Strong consistency is found between the theoretical calculations and n umerical results. Notably, with increasinga n dgradually increase. In practice, we prefer to use fewer PMLs for computational efficiency; however,the positive correlation betweend0and |R| is quite unfavorable for achieving this goal because using fewer PMLs corresponds to largerd0, resulting in worse absorption.

In practical applications, 10 PMLs are commonly used.Figure 12 shows a comparison among the reflections ofvxfor the PML+Dirichlet boundary, PML+UABC and PML+SABC for 10 PMLs. ForR0=7, the energy of reflection 1 dominates and focuses mainly atx= 670. ForR0=3, the energy of reflection 2 dominates and focuses mainly atx= 690. Thus, we use the maximum amplitude ofvxnearx= 690 to calculateandyielding 2.65 and 0.50, respectively, consistent with the results in Figure 11. The UABC amplifies the reflection considerably and the SABC reduces the reflection by only one-half.

Figure 11. Comparison between |R| derived from the semidiscrete plane wave solution and |R| calculated by numerical tests for the UABC and the SABC as a function of d0 with 30 PMLs.

Figure 12. Comparison among the reflections of vx for the PML+Dirichlet boundary, PML+UABC and PML+SABC with 10 PMLs. The reflections at approximately x = 670 and x = 690 are due to ε2 and ε1, respectively. Parts of the reflections overlap, indicating that using fewer PMLs leads to greater overlap.

Overall, the UABC amplifies reflection 2 unless a somewhat smalld0is used, which is not practical; the SABC can weaken reflection 2, but its absorption ability is relatively weak in our tests, especially when fewer PMLs are employed. Nevertheless, we find that simply increasingR0can significantly reduce the reflections.

Figure 13 compares the reflections among three combinations ofvxfor 30 PMLs withR0=4. In this case,reflection 1 dominates, and imposing the SABC is less meaningful. Moreover, for the PML+Dirichlet boundary withR0=3 (Figure 10), the amplitude of reflection 1 is 14.6, and that of reflection 2 is 113.3; for the PML+Dirichlet boundary withR0=4 (Figure 13), the amplitude of reflection 1 is 17.8, and that of reflection 2 is 10.2. ChangingR0from 3 to 4 (increasingd0from 32.2 to 39.2) causes the amplitude of reflection 2 to drop by 0.09 and that of reflection 1 to increase by 1.2. ForR0=3,although replacing the Dirichlet boundary with the SABC outside of the PML does not increase the amplitude of reflection 1, it makes the amplitude of reflection 2 drop by only 0.19. Thus, increasingR0from 3 to 4 is a better choice than replacing the Dirichlet boundary with the SABC, as the former not only improves the absorption of reflection 2 but also simplifies the implementation. The conclusion is the same when using 10 PMLs.

Figure 13. Comparison among the reflections of vx for the PML+Dirichlet boundary, PML+UABC and PML+SABC with R0 = 4; 30 PMLs are used.

Notably, although increasingR0from 3 to 4 can significantly weaken reflection 2, it slightly amplifies reflection 1. To clearly show the relations betweend0and the two reflections, we conduct a simple test. Figure 14 shows the maximum amplitudes of reflection 1 from the PML-interior interface and reflection 2 caused by the PML exterior boundary as a function ofd0. Asd0increases,reflection 1 only increases linearly; however, reflection 2 exponentially decreases. Generally, reflection 1 is weak,and increasingR0(d0) does not prominently amplify reflection 1; however, it can exponentially decrease the amplitude of reflection 2. Thus, increasingR0is a good strategy to mitigate reflections in practice.

Figure 14. Maximum amplitudes of reflection 1 from the PML-interior interface and reflection 2 caused by the PML exterior boundary as a function of d0 with 30 PMLs.

4.4. Absorption performance of the PML+UABC,PML+SABC, and PML+Dirichlet boundary with R0 =4 in 3D tests

In this subsection, we compare the absorption performance of the PML+UABC and PML+SABC by employing two 3D models and perform further examination by increasingR0from 3 to 4.

The medium has a compressional-wave velocityvP=3000m/s, a shear-wave velocityvS=2000m/s and a density ρ=1500kg/m3. The grid spacing is 100 m, the total number of grids is 400 × 200 × 400, and the size of the model is 39.9 km × 19.9 km × 39.9 km. Ten standard PMLs (α =0,β=1) with the quadratic damping profile are used.

The source time function is a Ricker wavelet with a center frequency of 2 Hz and a delay of 0.5 s. Spatial Gaussian smoothing is used for the source with a time interval of 0.02 s and a total of 300 time steps. Here, we use a sixth-order staggered FD operator to calculate the first spatial derivatives inside the inner (physical) domain and gradually lower the order at the boundary to the second order to implement the different exterior boundary conditions. We adopt the fourth-order Runge-Kutta scheme for the time-marching method.

For convenience, we define the reflection coefficients caused by changingR0in the numerical tests as:

4.4.1 Model 1 for small incident angles

An explosive source is applied at (20 km, 10 km, 20 km).Figure 15 shows the locations of the source and the receiver line for model 1. The receiver line is oriented in the same horizontal xoy plane as the source at (y= 18.5 km,z= 20 km) and is separated from the PML-interior interface by 5 grids. The reference solution is calculated by a large model with the same parameters.

Figure 15. Locations of the source and the receiver line for model 1. The red star represents the location of the source, the black dashed line represents the receiver line, and the red dotted line represents the PML-interior interface. The blue triangles represent the two receivers on the receiver line.Receiver 1 is at (20 km, 18.5 km, 20 km), and receiver 2 is at(29 km, 18.5 km, 20 km).

To avoid interference from incident waves, we calculate the reflections by taking the total wavefield minus the reference solution. Figure 16 shows the reflections ofvxalong the receiver line for the PML+ different exterior boundaries and differentR0. In Figure 16a,reflections 1 and 2 are both clearly visible. Reflection 1 weakens rapidly with increasing incident angle and appears only at nearly normal angles of incidence. ForR0=3, as shown in Figures 16b and 16c, the PML+UABC amplifies reflection 2, whereas the PML+SABC weakens reflection 2. Comparing Figure 16d with Figure 16a suggests that increasingR0from 3 to 4 does not notably increase the amplitude of reflection 1 but rapidly decreases that of reflection 2. As shown in Figures 16c and 16d, the PML+Dirichlet boundary withR0=4 achieves better absorption than the PML+SABC withR0=3 at small angles of incidence, and both achieve almost the same absorption at other incident angles.

Figure 17 shows the seismograms of the total wavefield and the reflections for the PML+different exterior boundaries and differentR0at receivers 1 and 2. These plots further confirm the conclusions drawn from Figure 16.Specifically, as shown in Figure 17c,｜andare 1.67, 0.47 and 0.11, respectively. In Figure 17d,a ndare 1.29, 0.18 and 0.13,respectively. The UABC amplifies the reflections at the two receivers. The SABC weakens the reflections but is severely influenced by the incident angle. Thus, increasingR0from 3 to 4 can yield the best and most stable absorption at small incident angles.

Figure 16. Reflections of vx along the receiver line for the PML+ different exterior boundaries and different R0. (a) PML+Dirichlet, R0 = 3; (b) PML+UABC, R0 = 3; (c) PML+SABC, R0 = 3; (d) PML+Dirichlet, R0 = 4. The same colorbar is used in(a)-(d).

Figure 17. Seismograms of the total wavefield and reflections of vx for the PML+different exterior boundaries and different R0 at receivers 1 and 2. (a) Total wavefield at receiver 1; (b) total wavefield at receiver 2; (c) reflection at receiver 1; (d) reflection at receiver 2.

Figure 18. |R3 UABC|, |R3 SABC|, and |R3 4| derived from the semidiscrete plane wave solution and |R| calculated by numerical tests as a function of the incident angle. The blue circles represent the absorption by increasing R0 from 3 to 4 for the PML+Dirichlet boundary calculated by numerical tests.

Overall, for small incident angles, as the incident angle increases, reflection 1 weakens rapidly, and increasingR0from 3 to 4 only slightly amplifies reflection 1. The UABC always amplifies reflection 2. The SABC can weaken reflection 2, but its absorption performance is quite limited. Conversely, increasingR0from 3 to 4 can significantly decrease the amplitude of reflection 2.

4.4.2 Model 2 for large incident angles

An explosive source is applied at (20 km, 17.9 km,20 km). Figure 19 shows the locations of the source and the receiver line for model 2. The receiver line is oriented in the same horizontal xoy plane of the source at (y= 18 km,z= 20 km) and is separated from the PML-interior interface by 10 grids. The reference solution is calculated by a large model with the same parameters.

Figure 20 shows the reflections ofvxalong the receiver line for the PML+different exterior boundaries and differentR0. Because the source is close to the PML, the incident angle rapidly increases as the distance increases.In this case, reflection 1 is negligible, and reflection 2 dominates. As shown in Figures 20a-d, the conclusion is the same as that for model 1.

Figure 19. Locations of the source and the receiver line of model 2. The red star represents the location of the source, the black dashed line represents the receiver line, and the red dotted line represents the PML-interior interface. The blue triangles represent the two receivers on the receiver line.Receiver 1 is at (22 km, 18 km, 20 km), and receiver 2 is at(29 km, 18 km, 20 km).

Figure 20. Reflections of vx along the receiver line for the PML+different exterior boundaries and different R0. (a) PML+Dirichlet, R0 = 3; (b) PML+UABC, R0 = 3; (c) PML+SABC, R0 = 3; (d) PML+Dirichlet, R0 = 4. The same colorbar is used in(a)-(d).

Figure 21. |R3 UABC|, |R3 SABC|, and |R3 4| derived from the semidiscrete plane wave solution and |R| calculated by numerical tests as a function of the incident angle. The blue circles represent the absorption by increasing R0 from 3 to 4 for the PML+Dirichlet boundary calculated by numerical tests.

Overall, for large incident angles, reflection 1 is almost negligible, and reflection 2 dominates. The UABC can only slightly weaken reflection 2 at a large incident angle and amplifies reflection 2 at other incident angles. Thus,the UABC is not recommended in any situation. In contrast, the SABC can weaken reflection 2, but its absorption performance is relatively weak, exhibiting poor absorption at both small and large incident angles and working well only near a certain incident angle. Moreover,the SABC requires complex procedures; thus, replacing the Dirichlet boundary with the SABC outside of the PML is not worthwhile. Conversely, increasingR0from 3 to 4 can easily weaken reflection 2 at small incident angles and achieves an absorption performance comparable to that of the SABC withR0=3 at large incident angles.

5. Conclusions

In this paper, we systematically compare the absorption performance of hybrid PMLs combining the standard PML with the first-order unstretched Higdon ABC and the first-order stretched Higdon ABC.

Unlike the previous rough studies based on the total wavefield or the total reflections of the PML+UABC or PML+SABC, we first use the stain algorithm to verify two main kinds of reflections in the discretized PML problem:reflections from the PML-interior interface and reflections caused by the exterior PML boundary. Next, we study the reflection coefficients of the UABC and SABC by removing the absorption of the PML and avoiding the interference of other possible reflections. We derive the reflection coefficients in continuous space, the results of which indicate that the UABC cannot absorb waves effectively and that the SABC works well. Next, to study the effects of the PML parameters in discrete space, we further derive the reflection coefficients in semidiscrete space, the results of which conform to the outcomes of our 1D and 3D numerical tests. With increasingR0(d0), the reflection coefficients of the UABC and SABC also increase. This property is unfavorable for achieving efficient absorption by using fewer PMLs in practice,which involves a largerd0.

We then quantify the absorption performance of both the UABC and the SABC in detail. With increasing incident angle, the reflection coefficients of the UABC and the SABC both first decrease and then increase. Notably,the UABC amplifies the reflection in most cases and weakens the reflection only slightly near large incident angles. Thus, the PML+UABC is not recommended in any situation. Conversely, the SABC always weakens the reflection; however, it works well only near a certain incident angle and achieves poor absorption at other incident angles. Furthermore, the SABC requires additional complex procedures; thus, replacing the Dirichlet boundary with the SABC outside of the PML is not worthwhile.

IncreasingR0(d0) can weaken reflections for a broad range of incident angles. This does not require any additional procedures; rather, it requires only a simple adjustment to the PML parameters. We increaseR0from 3 to 4 in our tests to facilitate a comparison and find that this adjustment results in better absorption than that of the PML+SABC. Thus, this scheme is recommended in practical applications; in fact,d0can be increased even larger for better absorption when the incident angle is large.

Our results show that the improvement in absorption achieved by combining the PML with either the SABC or the UABC is not better than that obtained by simply adjusting the damping profile of the PML; thus, the combination of the PML with the ABC is not recommended in practice.

Acknowledgments

This study was supported by the National Key R&D Program of China (No. 2018YFC1504204), National Natural Science Foundation of China (No. U1901602),Key Special Project for Introduced Talents Team of Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (No. GML2019ZD0203),Shenzhen Science and Technology Program (No.KQTD20170810111725321), and Shenzhen Key Laboratory of Deep Offshore Oil and Gas Exploration Technology (No. ZDSYS20190902093007855). The computational work was supported by the Center for Computational Science and Engineering of Southern University of Science and Technology.