Comparison of artificial absorbing boundaries for acoustic wave equation modelling

Key words: absorbing boundary condition, Higdon boundary, hybrid ABC, perfectly matched layer, sponge boundary.

Numerical modelling of seismic wavefields is important for understanding wave phenomena in complex media and is essential for full waveform inversion. Due to the restrictions on both memory requirement and computational cost, the model has to be limited in size and focused on the area of interest by introducing artificial boundaries. Therefore, an artificial boundary condition is needed to absorb the energy of the reflections from these artificial boundaries. Two main kinds of solutions have been proposed for this purpose: absorbing boundary conditions (ABCs) (e.g. Clayton and Engquist, 1977; Reynolds, 1978; Liao et al., 1984; Higdon, 1986; Higdon, 1991) and absorbing boundary layers (e.g. Cerjan et al., 1985; Kosloff and Kosloff, 1986; Compani-Tabrizi, 1986, Sochacki et al., 1987; Bérenger, 1994; Komatitsch and Martin, 2007; Liu and Sen, 2010, 2012).

The ABC splits the wave equation into two directions: inside and outside equations using the one-way wave equation method (Claerbout, 1985); then, only the outside equation is used on one or two layers outside of the interested area to avoid reflections inwards (as shown in Figure 1a). The ABC has a good performance when the incident waves are propagating within a certain angle range, especially when the incident waves are close to the direction normal to the boundary. When the incident angle is large (say 60° away from the direction normal to the boundary), the energy of artificial reflections would be boosted because the accuracy of both travel time and amplitude calculated with the one-way wave equation at wide angles is low.

Fig. 1.  Sketch map of the (a) absorbing boundary condition and (b) absorbing boundary layers. The square area in the middle is working area. The thick dashed-dot line in (a) represents the ABC applied outside the working area. The black solid lines in (b) represent the absorbing boundary layers applied outside the working area.

The ABC has been popular since the early 1970s, and some of the classical ABCs can be defined up to any desired order; however, the appearance of increasingly high-order derivatives in these ABCs renders them impractical beyond a certain order, typically 2 or 3. For example, the m-order Higdon boundary involves m-order derivatives both in space and time, and is thus very inconvenient for implementation when m is large (Bécache et al., 2010). Collino (1993) devised a new scheme for high-order ABCs by using special auxiliary variables. This scheme is based on a reformulation of the sequence of ABC that was proposed by Higdon (1986). In contrast to the original formulation of the Higdon conditions, this scheme does not involve any high derivatives beyond the second order by introducing special auxiliary variables. As a result, this scheme can be easily used up to any desired order m. Moreover, the computational cost only increases linearly with m (Givoli and Neta, 2003; Givoli, 2004). Corner compatibility conditions are derived for high-order radiation boundary conditions using auxiliary variable equations (Hagstrom and Warburton, 2004; Bécache et al., 2010).

The absorbing boundary layers use many layers to attenuate the artificial reflections gradually; thus, we can greatly reduce artificial reflections with an adequate number of layers (as shown in Figure 1b). There are three kinds of absorbing boundary layers: the sponge boundary (e.g. Cerjan et al., 1985; Compani-Tabrizi, 1986; Kosloff and Kosloff, 1986; Sochacki et al., 1987), the perfectly matched layer (PML) (e.g. Bérenger, 1994; Chew and Liu, 1996; Hastings et al., 1996; Collino and Tsogka, 2001; Marcinkovich and Olsen, 2003; Wang and Tang, 2003), and the hybrid absorbing boundary condition (Liu and Sen, 2010, 2012).

The sponge boundary (e.g. Cerjan et al., 1985; Compani-Tabrizi, 1986; Kosloff and Kosloff, 1986; Sochacki et al., 1987) is composed of dozens of layers outside the interested area. The sponge boundary introduced by Cerjan et al. (1985) avoids apparent inwards reflections by directly attenuating the wavefield with a gradually enhanced attenuation factor from inner to outer of the damping boundary belt. This method is simple in numerical implementation since there is no need to modify the wave equations. Sochacki et al. (1987) implemented the attenuation of seismic waves in the sponge boundary by introducing attenuation term directly to the wave equation. Zhou (1988) pointed out that a careful selection of the width of the damping boundary belt and attenuation factor is essential for the success of reducing the artificial reflections. However, no perfect rule for the selection has been available until now and one has to choose optimal parameters via numerical tests before it is ready for practical applications.

The PML (Bérenger, 1994) applies a completely new mechanism to avoid apparent artificial reflections. The PML introduces physical attenuations to the wave equation. The PML modifies the partial derivatives in the wave equation using complex coordinate stretching by introducing an imaginary part associated with an attenuation factor. Complex coordinate stretching is well known for viscoelastic media to understand the nature of the intrinsic attenuation for wave propagation; thus, the PML presents a nice rule in designing the optimal attenuation coefficients by tuning the attenuation factor. Rabinovich et al. (2010) compared high-order ABC with the PML in the frequency domain.

Although the traditional PML has been widely used in seismic wave simulations, it produces apparent artificial reflections for nearly grazing incident waves, low-frequency waves, and evanescent waves (e.g. Festa and Vilotte, 2005; Komatitsch and Martin, 2007; Drossaert and Giannopoulos, 2007a). Complex frequency-shifted PML (Kuzuoglu and Mittra, 1996) has a better absorbing performance in these cases (e.g. Festa et al., 2005; Festa and Vilotte, 2005; Drossaert and Giannopoulos, 2007a, 2007b; Komatitsch and Martin, 2007). In addition, the traditional PML adopts a non-physical splitting of wave equations, which has been proved to be weakly well posed (Abarbanel and Gottlieb, 1997). Thus, the convolutional PML is proposed to overcome this problem (e.g. Roden and Gedney, 2000; Komatitsch and Martin, 2007). Roden and Gedney (2000) further derived a non-split wavefield of convolutional complex frequency-shifted PML (CPML) that can be efficiently calculated with recursive convolution algorithm (Luebbers and Hunsberger, 1992). Kristek et al. (2009) discussed the split/unsplit, classical/convolutional and general/special PML formulations in detail, which can help the reader understand the classification of different PML formulations.

Recently, the PML with auxiliary differential equations (ADE-PML) has been proposed (e.g. Gedney and Zhao, 2010; Martin et al., 2010; Zhang and Shen, 2010). Both CPML and ADE-PML can be implemented for complex frequency-shifted operators, but the ADE form allows for extension to higher order time schemes and is also easier to implement for avoiding the convolutional calculation. Xie et al. (2014) gave a thorough review on the development from PML to ADE-PML.

Liu and Sen (2010, 2012) proposed a hybrid ABC. They split the model into three parts: the working area, the transition area and the ABC. The wavefields within the transition area are averaged by a linear weighting function between the wavefields generated by a two-way wave equation (i.e. acoustic wave equation) and one-way wave equation. The transition area smoothly absorbs the wavefields that propagate outside the working area and thus can decrease the boundary reflections. The hybrid ABC performs better than either the pure ABC or pure sponge boundary.

In the field of seismic exploration, artificial boundary conditions are becoming more and more important with an increasing demand on both the accuracy and computational efficiency of the migration and inversion, especially for large-scale 3D complex models. We need to further reduce computational cost and memory demands for the artificial boundary conditions. However, various artificial boundary conditions arise in history; thus, it is not so easy to be familiar with all of them, although this is necessary for further improvement towards much faster and smarter absorbing boundary conditions.

In this paper, we provide a thorough review of all typical absorbing boundary conditions and derive their equations in a uniform mathematical form. Then, we reveal the fundamental similarities and differences between similar methods through theoretical analyses. Next, we examine the performance of these boundary conditions via numerical experiments and qualitatively show their advantages and disadvantages. Finally, we summarise their applicable conditions and give some comprehensive suggestions on choosing different boundary conditions for practical applications.

Engquist and Majda (1977) discussed absorbing boundary conditions for a general class of differential equations based on pseudo-differential operators. Clayton and Engquist (1977) derived the absorbing boundary condition applicable to both acoustic and elastic wave equations by the approximation of a one-way wave equation. Engquist and Majda (1979) tried to solve the stability problem at the corners. Here we call this method the Clayton and Engquist boundary.

For a medium of the form D = {(x, z, t)|−a ≤ x ≤ a, −b ≤ z ≤ b}, a 2D scalar wave equation can be expressed as

We apply the Fourier transform on both temporal and spatial variables to Equation 1, and the dispersion relation can be obtained as follows

where k_x and k_z are the horizontal and vertical wavenumbers, respectively; ω is the circular frequency, p is the displacement, v is the velocity of the acoustic wave propagating in the media. According to the one-way wave equation theory (Claerbout, 1985), we can obtain the square-root operator along x-direction as follows

where ± represents a plane wave propagating along the positive or negative x-axis, respectively. The one-way wave equation method can separate the incoming and outgoing wavefields around the boundaries. For the positive x-axis at the right boundary, we have

which is the dispersion relation of the one-way wave equation that controls the wave propagations outwards. The one-way wave equation is solved iteratively along spatial directions; thus, only the wavefield on the previous one or two layers are needed to generate the wavefield in the current layer. In other words, wavefields beyond the boundary are not needed, which means no boundary reflection occurs at the outer layer (i.e. the right boundary).

Taylor-series expansion of the square-root operator on the right side of Equation 4 leads to unstable differencing schemes (Engquist and Majda, 1977); thus, Clayton and Engquist (1977) expanded this operator by Padé approximation. The first three orders are as follows (Clayton and Engquist, 1977):

and

The recurrence relation for the above approximations is

where a₁ = 1.

Equations 6 and 7 are widely used in designing one-way wave equation migration algorithms (Claerbout, 1985). The differential formats corresponding to Equations 5–7 are as follows (Clayton and Engquist, 1977)

and

where A1, A2, and A3 are right absorbing boundary conditions using the first, second and third order paraxial approximation, respectively. The accuracy analyses of these three approximations are shown in Figure 2, compared with the analytical solution. Obviously, high orders have better accuracy than lower orders do. But the error of the third order is still apparent for wide angles away from the normal direction to the boundary. On the other hand, it is a little difficult to solve the third order approximation shown in A3; in contrast, A1 is simple in implementation but its accuracy is too low. Therefore, A2 is a reasonable trade-off and is widely used. The differential equations for A2 on the left, upper and bottom boundaries are as follows (Clayton and Engquist, 1977):

Fig. 2.  Dispersion relations for the scalar wave equation. Curves A1, A2 and A3 are the dispersion relations of the first three orders using paraxial approximations to the scalar wave equation, respectively.

The Clayton and Engquist boundary needs the value of wavefields both along x-axis and z-axis; thus, the corners would act as point sources and cause instabilities (Engquist and Majda, 1979). Therefore, special treatments are needed at the corners for high-order Clayton and Engquist boundary conditions. One can replace A2 with A1 at the corners by assuming that the incident angle is 45° at the corners to avoid numerical instabilities (Clayton and Engquist, 1977; Engquist and Majda, 1979). For example, the differential formula for the lower right-hand corner becomes (Clayton and Engquist, 1977)

Similarly, the other three corners for the absorbing boundary conditions are expressed as

According to the plane wave solutions of the wave equation, plane waves spread to the right boundary can be written as

where θ is the incident angle, namely the included angle between the wavefront and the x-axis; k_x and k_z are the wavenumbers along x-axis and z-axis, respectively; and k is the wavenumber along the incident angle. Similarly, the reflected wave can be expressed as

where r is the reflection coefficient. The total wavefields around the boundary are

Taking Equation 17 into the right boundary conditions, we can obtain the reflection coefficient (Clayton and Engquist, 1977) as

where j = 1, 2, 3 represents the first, second and third order paraxial approximation, respectively; and θ is the incident angle. Figure 3 shows the absolute value of the reflection coefficient versus incident angle for different orders of the Clayton and Engquist boundaries.

Fig. 3.  Comparison of the absolute value of the reflection coefficient for the Clayton and Engquist, Reynolds, and Higdon boundaries.

Reynolds (1978) derived another absorbing boundary condition based on the one-way wave equation, which is known as the transparent boundary condition. Here we call it the Reynolds boundary. The wave Equation 1 can be rewritten as follows:

Denoting

Equation 19 becomes

where 1/v ∂p/∂t + L₁p = 0 corresponds to the right boundary, and 1/v ∂p/∂t−L₁p = 0 corresponds to the left boundary. Using Taylor-series expansion L₁ can be approximated as follows:

Thus, on the right boundary we have

Substituting

into Equation 23, we obtain

Equation 25 is the expression of the Reynolds boundary on the right side.

Considering the stability condition of the finite-difference method, we can replace the coefficient of ∂²p/∂z² in Equation 23 by s/(1+s), where s = vΔt/Δx, Δt is the temporal interval, and Δx is the spatial interval in x-direction; thus, Equation 23 becomes (Reynolds, 1978)

which can be rearranged as

Obviously, Equation 23 is the special case when s = vΔt/Δx = 1.

Similarly, the Reynolds boundaries on the left, upper, and bottom boundaries are as follows (Reynolds, 1978):

Taking the right boundary as an example, we substitute Equation 17 into Equation 26 and obtain

where θ is the incident angle, M = j(ωt−kx cos θ±kz sin θ), and N = j(ωt + kx cos θ ± kz sin θ). Thus, the reflection coefficient of the Reynolds boundary versus the incident angle is (Reynolds, 1978)

Figure 3 shows the absolute value of the reflection coefficient versus incident angle using different s for the Reynolds boundary.

Bayliss and Turkel (1980) proposed an absorbing boundary condition in the spherical coordinate system. Higdon (1986) further derived a progressive absorbing boundary equation in the Cartesian coordinate system (Higdon, 1986, 1987, 1990, 1991). We call this absorbing boundary condition the Higdon boundary for simplicity. Peng and Toksöz (1995) optimised the coefficients to improve the absorbing performance.

Plane waves spread to the right boundary can be written as Equation 15, which also satisfies the boundary condition

Equation 31 can be regarded as a compatibility condition for Equation 15, and it can annihilate waves moving at angles of incidence ±θ. In particular, the boundary condition ∂p/∂x + 1/v ∂p/∂t = 0 is compatible with outgoing waves moving at normal incidence, which is the first-order Clayton and Engquist boundary. Similarly, a linear combination of plane waves moving outward at angles ±θ₁,…, ±θ_m would exactly satisfy the high-order versions (Higdon, 1986)

where m is the order of the absorbing boundary condition, and θ_j is the angular parameter that can be chosen. Higdon (1986) demonstrated that the application of Equation 32 as the absorbing boundary can completely absorb the incident waves along the incident angle ±θ_j in theory.

Higdon (1991) pointed out that the angle θ_j can be chosen to take advantage of a priori information about the direction that along which waves are expected to approach the boundary. For example, if waves near normal incidence are of greatest importance, then we use θ_j = 0 for all j. If a wide range of incident angles presents, then it may be advisable to move θ_j away from 0° so as to distribute the roots of the reflection coefficient at a broad range of angles. In general, the optimal choice of θ_j is problem dependent. We use θ₁ = 0 and θ₂ = π/6 as the angle parameters for the second-order Higdon boundary in this paper.

Similarly, we can get the equations on the left, upper, and bottom boundaries as follows (Higdon, 1991):

When m = 1, 2, 3, Equation 32 can be written as follows (Higdon, 1991)

and

B1, B2 and B3 are the first, second and third order Higdon boundaries on the right side, respectively. Similar to the Clayton and Engquist boundary, it is complex to solve B3 but B1 is too simple; thus, B2 is preferred for practical applications. Substituting Equation 17 into Equation 32, we obtain the reflection coefficient of the Higdon boundary (Higdon, 1986)

where θ is the incident angle, and θ_j is a preferred incident angle that hopes to be perfectly absorbed. Figure 3 shows the absolute value of the reflection coefficient of the Higdon boundary versus angles for different orders.

For the right boundary of the second-order Higdon boundary, the differential operator form of Equation 32 when m = 2 is as follows (Higdon, 1987; Sun, 2003)

where I is the unit operator, which satisfies E_xI = E_x and AI = A; A and B are the control parameters for the finite-difference scheme. Differential operator expressions on space and time are as follows:

where i and j are indices of spatial grids along x- and z-directions, respectively; and k is the index of temporary grids.

Equation 38 includes three type of finite-difference schemes: the forward Euler when A = 0, B = 1, the backward Euler when A = 0, B = 0, and the box scheme when A = 0.5, B = 0.5. We use the forward Euler difference scheme in this paper. For more details of other cases, readers can refer to Higdon (1987).

Equation 25 is the Reynolds boundary on the right side, under the special case when s = vΔt/Δx = 1. Equation 10 is the second-order Clayton and Engquist boundary on the right side. According to Equation 1, we replace ∂²p/∂z² in Equation 10 by 1/v²∂²p/∂t² − ∂²p/∂x² and obtain

We see that Equation 40 is exactly the same as Equation 25; thus, the second-order Clayton and Engquist boundary is the special case when s = vΔt/Δx = 1 of the Reynolds boundary.

As shown in Figure 3, when the incident angle is less than 45°, the absolute value of the reflection coefficient of the second-order Clayton and Engquist boundary is less than that of the Reynolds boundary for s = 0.25 and s = 0.5, but the absolute value of the reflection coefficient of the Reynolds boundary is large when s = 0.75. From 45° to 60°, the absolute value of the reflection coefficient of the second-order Clayton and Engquist boundary is less than that of the Reynolds boundary when s = 0.25 but larger when s = 0.5 and s = 0.75. When the incident angle is greater than 60°, the absolute value of the reflection coefficient of the second-order Clayton and Engquist boundary is larger than those of the Reynolds boundaries.

At the four corner points, high-order (from the second-order) Clayton and Engquist boundaries need special treatments, which are derived by assuming that the incident angle is 45°. However, in actual simulation, especially for media with complex structures, we do not know the exact direction of the incident wave; thus, the Clayton and Engquist boundary would encounter strong reflections from corners when incident waves are far away from 45°. In contrast, the Reynolds boundary does not need any special treatment at corner points and is more straightforward in implementation.

From Equations 9–11, and Equations 34–36, we see that A1 and B1 are the same when cosθ₁ = 1; A2 and B2 are the same when cosθ₁ = cosθ₂ = 1; A3 and B3 are the same when cosθ₁ = cosθ₂ = cosθ₃ = 1. Therefore, the Clayton and Engquist boundary is the special case of the Higdon boundary when the incident angle used is equal to 0°.

As shown in Figure 3, the absolute value of the reflection coefficient of the first-order Higdon boundary is exactly the same as that of the first-order Clayton and Engquist boundary. The second-order Higdon boundary always has a smaller absolute value for the reflection coefficient than the second-order Clayton and Engquist boundary does. Similarly, the third order Higdon boundary always has smaller absolute value of the reflection coefficient than the third order Clayton and Engquist boundary does.

For the second-order boundary, the Reynolds boundary seems to have much smaller absolute value of the reflection coefficient, ~60° compared to both the Clayton and Engquist boundary and the Higdon boundary; however, its absolute value of the reflection coefficient is much bigger than the latter, between 20° to 40°. This means that the actual performance of the Reynolds boundary would not be as good as that of the other two methods.

Cerjan et al. (1985) proposed to eliminate the reflected wave by setting the damping boundary with multilayers outside the working area. An outward decaying exponential function can effectively reduce incident waves. In each time step of the wavefield extrapolation, the amplitude of the seismic wave on each grid within the absorption area is decayed by the Gauss function (Cerjan et al., 1985)

where λ is the attenuation coefficient, n is the grid index of the damping area, and n₀ is the layer number of the damping area (i.e. the index of the outer most layer). Cerjan et al. (1985) gave a set of empirical parameters n₀ = 20 and λ = 0.015 through numerical experiments. Bording (2004) pointed out that these parameters are not optimised and suggested a group of optimised parameters n₀ = 45 and λ = 0.0053.

Cerjan et al. (1985) showed that decaying the amplitudes of the waves by multiplying by an exponential function in a region surrounding the model can considerably reduce reflections for waves at any angle of incidence. Their method, however, acts on discrete numerical solutions rather than on the wave equation. Sochacki et al. (1987) implemented the attenuation of seismic waves in the sponge boundary by introducing an attenuation term directly to the wave equation. The wave equation is extended as (Sochacki et al., 1987)

where A(x, z) is the attenuation function. In the working area, A(x, z) = 0; in the damping area, Sochacki et al. (1987) suggested five kinds of damper as follows:

1. Linear damper (Sochacki et al., 1987)

where α = S/M_xN_z, x₁ = x₀ + Δx, z₁ = z₀ + Δz, M_x is the length of attenuation zone along the x-direction, N_z is the length of the attenuation zone along the z-direction, and S = ‖‖ A ‖‖_∞ = max{A(x, z)};

2. Exponent damper (Sochacki et al., 1987)

where β = lnS/[lnM_xN_z−lnΔxΔz], and α = (ΔxΔz)^−β;

3. Cubic damper (Sochacki et al., 1987)

where α = S/(M_xN_z)³;

4. Exponential damper (Sochacki et al., 1987)

where β = ln S/[(M_x−Δx)(N_z−Δz)], and α = exp(−βΔxΔz);

5. Gaussian damper (Sochacki et al., 1987)

where β = {(Δx+Δz)(M_x + N_z)/[(M_x + N_z)−(Δx+Δz)]}ln S, and α = exp[−β/(Δx + Δz)].

Compared with the Cerjan sponge boundary, the Sochacki sponge boundary has the following advantages: first, the damping term has a physical meaning and can be used to study the effects of friction and other dissipative phenomena on acoustic and elastic waves; second, it is easier to determine the reflection coefficient; third, since the damping term is part of the wave equation, a variety of numerical techniques can be used and theoretical results are more easily obtained.

Bérenger (1994) introduced a completely new boundary condition called the perfectly matched layer (PML). The PML has a remarkable property of having a zero reflection coefficient for all incident angles and all frequencies for the continuous wave equation. The PML shows great superiority over the classical boundary conditions mentioned above and is widely used in acoustic wave simulations (e.g. Abarbanel and Gottlieb, 1997; Liu and Tao, 1997; Qi and Geers, 1998; Liu, 1999; Katsibas and Antonopoulos, 2002; Diaz and Joly, 2006; Bermúdez et al., 2007) and elastic wave simulations (e.g. Chew and Liu, 1996; Hastings et al., 1996; Collino and Tsogka, 2001; Festa and Nielsen, 2003; Komatitsch and Tromp, 2003; Basu and Chopra, 2004; Festa et al., 2005, Festa and Vilotte, 2005; Appelö and Kreiss, 2006; Komatitsch and Martin, 2007; Yang et al., 2007; Lan et al., 2013). The PML has been further extended to other methods, such as the pseudo-spectral method (Liu, 1998), the finite element method (Collino and Tsogka, 2001), the spectral element method (Festa and Vilotte, 2005), and the grid method (Xu and Zhang, 2008).

The 2D scalar wave equation shown in Equation 1 can be rewritten as (Liu and Tao, 1997; Hustedt et al., 2004)

After transforming into the frequency domain by the Fourier transform, Equation 48 becomes

where P_x, P_z, P, , and are the Fourier transforms of P_x, P_z, p, A_x and A_z, respectively.

The essence of the PML is to replace propagating (oscillating) waves with exponentially decaying waves by analytically continuing the wave equation into complex coordinates (e.g. Bérenger, 1994; Collino and Tsogka, 2001; Nataf, 2013) using

where d_x (x) is the attenuation coefficients in x-direction, which can be obtained by (Collino and Tsogka, 2001)

where v_max is the maximum acoustic wave velocity, δ is the thickness of the PML area, and R is the theoretical reflection coefficient. We use R = 10⁻⁵ when δ = 10Δx, and R = 10⁻⁷ when δ = 20Δx.

Substituting Equation 50 into Equation 49, we obtain

After rearranging Equation 52, we have

Finally, the acoustic wave equation with PML can be expressed as (Liu and Tao, 1997)

The wave equation is extended to complex coordinates by (Komatitsch and Martin, 2007)

where α_x = α_max(1−x/δ), α_max = πf₀, f₀ is the dominant frequency of the source, and χ_x = 1 + (χ_max −1) (x/δ)². Generally χ_x = 1 is good enough for most applications (Komatitsch and Martin, 2007). Thus, we can obtain

On the assumption that is the inverse Fourier transform of 1/s_x about ω, we have

Denoting

we obtain

In the time domain, we have

where and ∂_x = ∂/∂x. At the nth time step, the convolutional term in Equation 60 can be written as

The integral in Equation 61 is written as a discrete form:

Using Equations 61 and 62, we obtain

Let b_x = exp[ −(d_x/χ_x + α_x)mdt], a_x = (b_x −1)d_x/[χ_x(d_x + χ_x α_x)]. Using the recursive convolution algorithm deduced by Luebbers and Hunsberger (1992), we have

thus,

Pasalic and McGarry (2010) gave a detailed derivation of CPML for the second-order acoustic wave equation. Introducing new auxiliary variables φ_x and φ_z, we can rewrite the second-order spatial derivatives in terms of the first-order derivatives as (Pasalic and McGarry, 2010)

where

The CPML implementation for the acoustic wave equation is expressed as (Pasalic and McGarry, 2010)

Liu and Sen (2010, 2012) proposed an efficient hybrid ABC to absorb boundary reflections. The transition area is added between the working area and the model boundary. Figure 4 shows the sketch map of the hybrid ABC. The wavefields in the transition area (grey area in Figure 4) are weighted between two-way wavefields and one-way wavefields, where the one-way wavefields are calculated by the second-order Higdon boundary. A linear weighting function is used to sum the wavefields (Liu and Sen, 2010, 2012).

Fig. 4.  Sketch map of the hybrid absorbing boundary condition (ABC). The grey zone is the transition area.

where p^two are the two-way wavefields, p^one are the one-way wavefields, w_Ai = (i−1)/N, and N is the number of layers used in the hybrid ABC.

We perform numerical experiments on a homogeneous medium. The wave velocity of a square model is 2500 m/s. The spatial grid interval is 5 m, and the grid number is 301 × 301. The source is a Ricker wavelet with a dominant frequency of 20 Hz. We compare the boundary reflections using various boundary conditions: the second-order Clayton and Engquist boundary; the second-order Higdon boundary; the Reynolds boundary; the Cerjan sponge boundary (with 20 absorbing layers) (using parameters given by Cerjan et al., 1985); the Sochacki sponge boundary (20 layers) (using the linear damper given by Sochacki et al., 1987); the hybrid ABC (10 layers and 20 layers); the PML (10 layers and 20 layers); and the CPML (10 layers and 20 layers).

Artificial absorbing boundary conditions are used on all boundaries of the model. We perform 12 sets of experiments to compare the reflections of several absorbing boundary conditions. We use a very large model to simulate the theoretical wavefield to avoid boundary reflections, which can be regarded as a reference to check the performance of the artificial absorbing boundary conditions.

Figures 5 and 6 show the snapshots obtained by different boundary conditions. Figures 7 and 8 show the waveforms obtained by different boundary conditions. Figure 9 shows the wavefield records obtained on the ground surface (i.e. 0 m in depth). The maximum value of the wavefield records is 0.1294. Among the three ABC, the absorption performance of the Clayton and Engquist boundary is better than that of the Reynolds boundary, and the absorption performance of the Higdon boundary is the best. The hybrid ABC (Liu and Sen, 2010, 2012) is shown to be superior to the Higdon boundary. Compared with the PML and CPML, however, absorption performance of the hybrid ABC is relatively poor. In terms of the computational cost, at least, the Higdon boundary is superior to the absorbing layers (i.e. hybrid ABC, PML, CPML or sponge boundary) since only one or two layers are involved in the computations. Therefore, the Higdon boundary would be a good candidate when the demand of absorbing performance is not so serious but the computational cost is heavy.

Fig. 5.  Snapshots of the wavefield using different boundary conditions at 450 ms: (a) theoretical wavefield; (b) the Clayton and Engquist boundary; (c) the 20-layer Cerjan boundary (Cerjan-20); (d) the 10-layer hybrid absorbing boundary condition (ABC) (Hybrid-10); (e) the 10-layer perfectly matched layer (PML-10); (f) the 10-layer convolutional complex frequency-shifted perfectly matched layer (CPML-10); (g) the Reynolds boundary; (h) the Higdon boundary; (i) the 20-layer Sochachki boundary (Sochachki-20); (j) the 20-layer hybrid ABC (Hybrid-20); (k) the 20-layer PML (PML-20); and (l) the 20-layer CPML (CPML-20).

Fig. 6.  Snapshots of the wavefield using different boundary conditions at 650 ms: (a) theoretical wavefield; (b) the Clayton and Engquist boundary; (c) the 20-layer Cerjan boundary (Cerjan-20); (d) the 10-layer hybrid absorbing boundary condition (ABC) (Hybrid-10); (e) the 10-layer perfectly matched layer (PML-10); (f) the 10-layer convolutional complex frequency-shifted perfectly matched layer (CPML-10); (g) the Reynolds boundary; (h) the Higdon boundary; (i) the 20-layer Sochachki boundary (Sochachki-20); (j) the 20-layer hybrid ABC (Hybrid-20); (k) the 20-layer PML (PML-20); and (l) the 20-layer CPML (CPML-20).

Fig. 7.  Waveforms of boundary reflections using different absorbing boundary conditions at 650 ms. This is a horizontal slice in the snapshot crossing the source location.

Fig. 8.  Waveforms of boundary reflections using different absorbing boundary conditions at 650 ms. This is a horizontal slice in the snapshot crossing the source location. The enlarged portion of the black frame is shown in the figure.

Fig. 9.  Shot records of a point source with different absorbing boundary conditions: (a) theoretical wavefield; (b) the Clayton and Engquist boundary; (c) the 20-layer Cerjan boundary (Cerjan-20); (d) the 10-layer hybrid absorbing boundary condition (ABC) (Hybrid-10); (e) the 10-layer perfectly matched layer (PML-10); (f) the 10-layer convolutional complex frequency-shifted perfectly matched layer (CPML-10); (g) the Reynolds boundary; (h) the Higdon boundary; (i) the 20-layer Sochachki boundary (Sochachki-20); (j) the 20-layer hybrid ABC (Hybrid-20); (k) the 20-layer PML (PML-20); and (l) the 20-layer CPML (CPML-20).

The absorption performance of the sponge boundary is not as good as that of the PML method and even poorer than that of the Higdon boundary. Of course, the absorption performance of the sponge boundary may be improved after some careful tests and selections on the parameters (e.g. tuning the thickness of the damping belt and the attenuation coefficients). Nevertheless, our numerical results shown here reveal that the sponge boundary would have potential risks, and further research is still needed for reliable or universally optimal parameters.

As shown in Figure 8, the 20-layer PML (or CPML) performs better than the 10-layer one does, since the former has almost no visible boundary reflections. In addition, the 10-layer PML has smaller artificial reflections than the 20-layer hybrid ABC does. This indicates that the PML is superior to the hybrid ABC with the same number of absorbing layers.

In practical applications, our model is sometimes very large in scale, which would lead to a nearly grazing incidence when the wavefields are propagating far away from the source. We further examine various absorbing boundary conditions for a large-scale long model. The grid number is 601 × 61. The seismic source is located at 2500 m along x-direction and 50 m in depth. The other parameters are the same as the square model. Seven boundary conditions are compared: the second-order Higdon boundary, the hybrid ABC (with 10 and 20 absorbing layers), the PML (10 and 20 layers) and the CPML (10 and 20 layers). We use a much larger scale model to generate an artificial-reflection free case as the reference for checking the absorbing performance.

Figure 10 shows the snapshots of the wavefield at three different time steps. Figure 11 shows the differences between absorbing boundary conditions and an artificial-reflection free reference, which can be regarded as the error or artificial reflections introduced by the absorbing boundary condition. Figure 12 shows the waveforms recorded at a depth of 0 m at different time steps. The maximum value of the wavefield records is 0.1127. Obviously, the wavefront is moving away from the source position with increasing time, and the incident angle gradually increases; finally, the wavefront grazes the boundary at large distances, as shown in Figures 10 and 11. The Higdon boundary exhibits significant artificial reflections at grazing areas, as shown in Figure 11. Nevertheless, the hybrid ABC has a better absorbing performance than the Higdon boundary. However, the hybrid ABC is still less efficient for absorbing the nearly grazing incidence since it is based on a one-way wave equation.

Fig. 10.  Snapshots of a point source at 2500 m along x-direction and 0 m depth in the narrow slice model with different absorbing boundary conditions at 300 ms, 600 ms and 900 ms: (a) theoretical wavefield; (b) Higdon boundary; (c) Hybrid-10; (d) PML-10; (e) CPML-10; (f) Hybrid-20; (g) PML-20; and (h) CPML-20.

Fig. 11.  Snapshots of the wavefield difference of a point source at 2500 m along x-direction and 0 m depth in the narrow slice model with different absorbing boundary conditions at 300 ms, 600 ms and 900 ms: (a) theoretical difference; (b) wavefield difference between Higdon and theoretical; (c) wavefield difference between Hybrid-10 and theoretical; (d) wavefield difference between PML-10 and theoretical; (e) wavefield difference between CPML-10 and theoretical; (f) wavefield difference between Hybrid-20 and theoretical; (g) wavefield difference between PML-20 and theoretical; and (h) wavefield difference between CPML-20 and theoretical.

Fig. 12.  Waveforms using different boundary conditions along x-direction at 0 m depth: (a) a slice from 1750 m to 2375 m at 300 ms; (b) a slice from 250 m to 875 m at 900 ms.

In contrast, the PML performs better than the hybrid ABC, and the CPML performs better than the PML does. The 20-layer PML performs better than the 10-layer CPML when the grazing effects are not so obvious at short traveltimes; in contrast, when the grazing effects are strong at long travel times, the CPML shows it to be necessary to absorb these grazing waves, as shown in Figures 13 and 14. Figure 15 shows the single-trace records using different boundary conditions. For the convenience of comparing the absorption performances of different boundary conditions, we calculate the numerical reflection coefficient using

Fig. 13.  Shot records of a point source with different absorbing boundary conditions: (a) theoretical wavefield; (b) Higdon boundary; (c) Hybrid-10; (d) PML-10; (e) CPML-10; (f) Hybrid-20; (g) PML-20; and (h) CPML-20.

Fig. 14.  Wavefield difference of shot records with different absorbing boundary conditions: (a) theoretical difference; (b) wavefield difference between Higdon and theoretical; (c) wavefield difference between Hybrid-10 and theoretical; (d) wavefield difference between PML-10 and theoretical; (e) wavefield difference between CPML-10 and theoretical; (f) wavefield difference between Hybrid-20 and theoretical; (g) wavefield difference between PML-20 and theoretical; and (h) wavefield difference between CPML-20 and theoretical.

Fig. 15.  Single-trace records using different boundary conditions: (a) at 2500 m; (b) at 500 m.

where p_ref is the theoretical wavefield without artificial reflections, and the p_mod is the wavefield with different artificial boundary conditions. Figure 16a shows the numerical reflection coefficient for different boundary conditions at different incident angles, and Figure 16b shows the decibel (dB) values of numerical reflection coefficient calculated by

Fig. 16.  Comparison of the reflection coefficient computed by numerical simulations: (a) absolute values of numerical reflection coefficient; (b) the dB values of numerical reflection coefficient.

Figure 16b shows that the 10-layer CPML has a better absorbing performance compared with the 10-layer PML for the nearly grazing incident wave from 63° to 78°. However, the 20-layer CPML only has slight advantage over the 20-layer PML from 80° to 84°. With the same number of absorbing layers, the CPML has better absorbing performance than the PML does; however, we can still use more layers for the PML instead of using a thin CPML. In addition, we see that the 10-layer PML is better than the 20-layer hybrid ABC.

Therefore, the absorption performance from good to poor in sequence is: CPML-20, PML-20, CPML-10, PML-10, Hybrid-20, Hybrid-10, and Higdon.

We test the modified Marmousi model (Figure 17), whose grid spacing is 5 m, and the grid number is 737 × 751. The upper boundary of the model is a free surface, and the other three edges are absorbing boundaries. The source is a Ricker wavelet and the dominant frequency is of 10 Hz. The source is added on the free surface and 500 m away from the left-upper corner of the model. Three kinds of boundary conditions are compared numerically: the Higdon boundary, the hybrid ABC (10 and 20 layers), the PML (10 and 20 layers), and the CPML (10 and 20 layers). There is no theoretical wavefield available as a reference for the Marmousi model; thus, we use the 50-layer CPML as a reference instead. Two groups of receivers are located along the upper and bottom boundary, respectively. The interval between two adjacent receivers is 10 m.

Fig. 17.  Modified Marmousi velocity model.

Figures 18 and 19 show the wavefield records obtained by the upper and bottom receivers, respectively. The absorption performance of the Higdon boundary is poor compared with the other methods, especially for deeper reflections. This would be harmful for the full waveform inversion. Whereas, we can see that the waveforms obtained using the Higdon boundary and hybrid ABC shown in Figure 18 are fairly close to those obtained using the PML or CPML. In fact, the artificial reflections may not be so serious to the migration results as to the full waveform inversion since the imaging result is not so sensitive to the amplitude. Thus, the Higdon boundary is still a good candidate for migration velocity analyses based on reverse time migration, due to its simplicity and very low computational cost as well as memory demand.

Fig. 18.  Shot records and wavefield difference of a point source on the surface of Marmousi model. (a-1) Higdon boundary; (b-1) Hybrid-10; (c-1) PML-10; (d-1) CPML-10; (e-1) Hybrid-20; (f-1) PML-20; (g-1) CPML-20; (a-2) wavefield difference between Higdon and CPML-50; (b-2) wavefield difference between Hybrid-10 and CPML-50; (c-2) wavefield difference between PML-10 and CPML-50; (d-2) wavefield difference between CPML-10 and CPML-50; (e-2) wavefield difference between Hybrid-20 and CPML-50; (f-2) wavefield difference between PML-20 and CPML-50; and (g-2) wavefield difference between CPML-20 and CPML-50.

Fig. 19.  Shot records and wavefield difference of a point source at the bottom of Marmousi model. (a-1) Higdon boundary; (b-1) Hybrid-10; (c-1) PML-10; (d-1) CPML-10; (e-1) Hybrid-20; (f-1) PML-20; (g-1) CPML-20; (a-2) wavefield difference between Higdon and CPML-50; (b-2) wavefield difference between Hybrid-10 and CPML-50; (c-2) wavefield difference between PML-10 and CPML-50; (d-2) wavefield difference between CPML-10 and CPML-50; (e-2) wavefield difference between Hybrid-20 and CPML-50; (f-2) wavefield difference between PML-20 and CPML-50; (g-2) wavefield difference between CPML-20 and CPML-50.

The absorption performance of the 20-layer CPML is the best among the three methods listed; however, its computational cost and memory demand are bigger than that of the 20-layer PML. In addition, it only shows negligible superiority over the 20-layer PML for most cases. Therefore, we can use 20-layer PML instead in most of the middle-sized models (i.e. the usual cases), except that the source is close to the left or right boundaries, which may encounter severe grazing effects.

In the above three numerical experiments, we use θ₁ = 0 and θ₂ = π/6 as the angle parameters for the second-order Higdon boundary. The main consideration is that the accuracy of the approximation of a one-way wave equation is generally no more than π/3 (60°) (Claerbout, 1985; Zhang et al., 2010); thus, we should not choose large angles. In complex media, we cannot know the exact incidence angle at the boundary; thus, we only guarantee that incident angles smaller than 60° can be handled to reduce the boundary reflections. We use θ₁ = 0 just as the Clayton and Engquist boundary does. Then, we select θ₂ = π/6 (30°), which is between 0° and 60°. Numerical simulations show that the absorption effect of the second-order Higdon boundary is always better than that of the second-order Clayton and Engquist boundary for a wider angle range. Higher order Higdon boundary conditions have more controllable angle parameters thus are more powerful in handling wider incident angles. However, higher order Higdon boundary conditions would have limited improvement due to the phase and the amplitude errors that are associated with the expansion of one-way wave equation (Claerbout, 1985).

Among the ABC, there is another method not mentioned in the text, which is called Liao’s ABC (Liao et al., 1984; Liao, 1996). Liao’s ABC extrapolates wave equations in the time and space domain based on the Newton’s backward differential. This method is seldom reported in seismic exploration but is widely used in electromagnetic simulations to absorb the boundary reflections (e.g. Wagner and Chew, 1995; Wang and Tang, 2010; Zhang and Yu, 2012). Wang and Tang (2010) point out that Liao’s ABC is closely related to Higdon’s space–time extrapolation method and has a nice feature of simplicity and clarity, which is of particular interest for high-order boundary treatments. To avoid confusion, we classify Liao’s method into the Higdon boundary condition since the Higdon boundary condition is already well known in the community, although Liao’s method (1984) was proposed two years earlier than the Higdon boundary (1986).

The artificial absorbing boundary condition plays an important role in the numerical simulation of seismic wavefields. Two main kinds of methods have been proposed for this purpose: ABCs and the absorption boundary layers (i.e. the sponge boundary, the hybrid ABC, and the PML/CPML). We derive similar absorbing boundary conditions using unified forms first, and then compare their absorbing performances by theoretical analyses and numerical experiments.

The Clayton and Engquist boundary has difficulties in dealing with the corner points, which may result in strong artificial reflections. The Reynolds boundary deals with corner points with ease, but has poor absorbing performance, since its tuning parameter would degrade the absorbing performance for the low angle incident waves when trying to improve the absorbing performance for other incident angles. The Clayton and Engquist boundary is basically a special case of the Higdon boundary; whereas, the Higdon boundary does not have problems at corner points and has the best absorbing performance among these three ABCs. Therefore, we suggest using the Higdon boundary at any time if one hopes to use the absorbing boundary condition, when the computational cost and memory demand are heavy but the demand for the absorbing performance is not so serious. For example, the migration is not so sensitive to the weak amplitude distortion caused by artificial reflections. In contrast, the full wave form inversion would have problems in the presence of weak amplitude distortion caused by artefacts.

The sponge boundary still has obstacles on determining the optimal parameters before a series of numerical tests. In addition, the sponge boundary usually has to use more layers (~30–50) to guarantee the absorbing performance; however, the PML with 10 or 20 layers would show great superiority over the sponge boundary both in computational cost and memory demand. Therefore, the sponge boundary is not suggested all the time for practical applications and one can use 10- or 20-layer PML instead.

The PML/CPML is not convenient to implement for the second-order formulation of acoustic and elastic wave equations. The hybrid ABC is the only method that can have comparable performance with the PML/CPML for the second-order wave equation. The hybrid ABC is easier in implementation than the PML and CPML, and it has a much better absorbing performance than the Higdon boundary. However, the hybrid ABC is still less efficient for absorbing nearly grazing waves since it is based on a one-way wave equation.

The 10-layer PML is accurate enough when the nearly grazing wave is not serious. The 20-layer PML is suggested for most practical applications for general size models, even in the presence of strong nearly grazing waves. The main advantage of PML is that it favours natural implementation of high-order temporal discretisation if necessary, while CPML does not. To implement CPML using high-order temporal discretisation, it is necessary to introduce auxiliary variables and equations, which will introduce extra costs in both time and memory. The 20-layer CPML is only necessary when the source is located fairly close to the boundary for large-scale models, where the nearly grazing wave is severe.