PML Absorbing Boundary Condition for Seismic Numerical Modeling by Convolutional Differentiator in Fluid-Saturated Porous Media

INTRODUCTION

Elastic wave propagation in a fluid-saturated po-rous medium is governed by Biot's theory (Biot, 1956a, b). In a homogeneous medium, analytical solu-tions for Biot's equations can be obtained easily, but such solutions are, in general, impossible for an arbi-trary heterogeneous medium. Elastic wave propaga-tion in complex fluid-saturated porous media is of great importance in solid geophysics such as seismol-ogy and other branches of applied geophysics such as petroleum engineering. Numerical methods are usually used to solve these problems. Examples of nu-merical methods include the finite-difference, finite-element, pseudospectral methods and convolutional methods. In this article, a Forsyte polynomial convo-lutional differentiator method is to be used.

In numerical simulation of seismic wave propa-gation in an unbounded medium, the introduction of artificial boundaries will introduce spurious reflec-tions that will affect the accuracy of numerical solu-tions. Although the problem can be overcome by in-creasing the size of the model, it is not always feasible because of the large amount of computer resource needed for a long-time solution. It is then highly de-sirable to eliminate these reflections by using absorb-ing boundary condition when modeling an unbounded medium. Obtaining an effective and stable absorbing boundary condition is always one of the most impor-tant tasks in numerical modeling of seismic wave propagation. Smith (1974) proposed an absorbing boundary condition for finite-difference and finite element methods. In his method, the Dirichlet and Neumann conditions are used alternatively, and the solutions from these two conditions are superimposed. Although it is easy to implement, this method greatly increases the computation time. A widely used ab-sorbing boundary condition in seismic modeling pro-posed by Clayton and Engquist (1977) is a one-way wave equation based on the paraxial approximations of the acoustic or elastic equations. Similar ap-proaches were proposed by several authors. Although effective for small incidence angles, these absorbing boundary conditions degrade for large angles of inci-dence. Another approach is to add a spatial filter or damping taper to the boundaries (Kosloff and Kosloff, 1986; Cerjan et al., 1985). In this so-called sponge absorber method, the transition zone from the inner region to the outer boundary should be thick and smooth. Liao and Wong (1984) developed an absorb-ing boundary condition based on the principle of plane wave and interpolation.

Berenger (1994) proposed a highly effective per-fectly matched layer (PML) as an absorbing boundary condition for electromagnetic waves. It has been widely used for finite-difference and finite-element methods (Liu, 1997; Chew and Weedon, 1994) since then. Chew and Liu (1996) first proved that such a PML also exists for elastic waves in spite of the cou-pling of S- and P-waves at an elastic interface. The PML has zero reflection to the regular elastic medium, although a small reflection can result from discretiza-tion in the PML scheme. The PML absorbing bound-ary condition has also been extended to model acous-tic waves and electromagnetic waves in lossy media (Liu, 1997; Liu and Tao, 1997) as well as electro-magnetic and elastic waves in cylindrical and spheri-cal coordinates (He and Liu, 1999; Liu, 1999).

These previous absorbing boundary conditions, however, have been developed for electromagnetic, acoustic, and elastic waves in solids. Little attention has been paid to elastic waves in fluid-saturated po-rous media. In this work, we extended the PML to truncating unbounded fluid-saturated porous media for numerical solutions using a Forsyte polynomial con-volutional differentiator method proposed by Li and Li (2008). We adopt the 2D homogeneous and het-erogeneous fluid-saturated porous media as examples to implement the numerical simulation by using the convolutional differentiator method with PML ab-sorbing boundary condition. Numerical results dem-onstrate that the PML absorbing boundary condition can work well with our method.

BASIC FORMULA

Biot's theory for elastic waves in porous media was established on a macroscopic level (Biot, 1956a, b). The following assumptions are necessary to the theory: (1) seismic wavelength is large in comparison to the pore size, (2) the deformations of the solid framework are small, (3) the fluid phase is continuous so that pores are connected, (4) the solid matrix is elastic, (5) the medium is statistically isotropic, and (6) gravity forces are neglected. Equations of motion in fluid-saturated porous medium are described as

where

By means of the parameters ρ₁₁, ρ₁₂ and ρ₂₂, we can rewrite equations (1) and (2) as

In fact, equations (4) and (5) represent the basic formula that the solid and fluid displacement compo-nent satisfies.

By discretization and transformation, the first-order velocity-stress wave equations can be given from the constitute equations in Biot media as follows.

The velocity component of the solid phase is

and the velocity component of the fluid phase is

while the stress components of the solid phase and fluid phase are

and

respectively. In the above equations, the parameters describing the physical property of the media are as follows: u_3×1 is average displacement vector of the solid phase and U_3×1 is average displacement vector of the fluid phase; b_i is dissipation coefficient and D_6×6 is elastic coefficient matrix in constitutive relations; P is normal stress of fluid material; v_i is particle velocity of solid material in i direction and V_i is particle velocity of fluid material in i direction; σ_ii is normal stress of solid material, while τ_ij(i≠j) is shear stress of solid material; and ρ₁₁ is density of the solid material, ρ₂₂ is density of the fluid material, and ρ₁₂ is coupling den-sity parameter. Q is the coupling parameter relating to bulk variation between the solid phase and the fluid phase, and R represents the elastic parameter relating to the fluid phase, while ε is the bulk strain of the fluid phase.

EQUATIONS SPLITTING IN TIME DOMAIN

The PML boundary conditions is implemented by setting a number of special grids surrounding the computation region to cause seismic waves to attenuate rapidly with increasing of the spatial grid until the wave could be absorbed completely.

At the time when we simulate seismic wave propagation by using the method put forward by Li and Li (2008), we should construct the different discrete schemes of the wave equations inside the model and in the PML region; the interface between the modeling region and the PML region will be treated as the inner interface.

Next, we will introduce step-by-step the PML boundary conditions for the 2D wave equations in two-phase media (no consideration of attenuation function).

First, we should separate the wavefield according to the x component and the z component; the velocity components are expressed as v_x = v_x^x + v_x^z, v_z = v_z^x + v_z^z, V_x = V_x^x + V_x^z, and V_z = V_z^x + V_z^z, while the stress components are expressed as σ_xx = σ_xx^x + σ_xx^z, σ_zz = σ_zz^x + σ_zz^z, τ_xz = τ_xz^x + τ_xz^z, and P = P^x + P^z.

In the above expressions, the superscript x de-notes that the wave propagates in x direction and the superscript z are the same.

Second, the wave equations are separated into the following parts:

The equations that the velocity of the solid phase satisfied are

and the equations that the velocity of the fluid phase satisfied are

Accordingly, the equations that the stress of the solid phase satisfied are

and the equations that the stress of the fluid phase satisfied are

Third, we will construct the discrete forms of the above-mentioned equations. The detailed discrete forms are omitted here for simplicity.

IMPLEMENTATION OF THE FPCD ALGORITHM

In the last section, we obtained 16 time-domain equations for 16 split field variables v_x^x, v_x^z, v_z^x, v_z^z, V_x^x, V_x^z, V_z^x, V_z^z, σ_xx^x, σ_xx^z, σ_zz^x, σ_zz^z, τ_xz^x, τ_xz^z, P^x, and P^z, which already incorporate the PML boundary condition. In numerical implementation, the time and space are discretized by time step Δt and node spacing Δx and Δz. The unbounded medium is truncated into a finite computational domain with a total of N_x×N_z grid nodes. Then, the FPCD method proposed by Li and Li (2008) is applied to each equation to obtain the updated variable at next step. In this CPFD scheme, the spatial and temporal derivatives are obtained by CPFD method and staggered-grid finite-difference method separately.

Specially, the PML layer is an added grid layer outside the computation region and a nonphysically virtual region. The PML layer is divided into four side boundaries and four corners, and the attenuation factors in different parts are different during the computation time.

The absorption of outgoing wave is achieved by the PML region, which consists of several cells of PML materials with a quadratically or linearly tapered profile to increase the attenuation toward the outer boundary. In this article, 20 cells of PML with a quadratic profile for d(x) and d(z) are used. Therefore, in the above-mentioned equations, the factors d(x) and d(z) have the following forms

where f₀ is the dominant frequency of the source, L_PML is the thickness of the PML cells, and lx_i is the distance from the interface between the interior region and PML region. After testing several values for a₀, we found that 1.79 has the best absorption of outgoing waves. The outer boundary conditions for the PML are chosen as a hard boundary or alternatively a soft boundary.

NUMERICAL EXAMPLES

We have implemented the Forsyte polynomial convolutional differentiator algorithm with the PMLs as an absorbing boundary condition for Biot's equa-tions in 2D fluid-saturated media. Unlike the continu-ous case, a small reflection will occur at the PML in-terface due to the discretization and truncation. To minimize the reflection from the PML region, we use a quadratic profile with 20 cells of PMLs at the com-putational edge. In the following examples, a line, compressional source is used to excite the seismic wave. The source time function is the Ricker wavelet

where f₀ is the dominant frequency, and t₀ is the cen-tral time of the wavelet.

A Homogeneous Model

In this experiment, a homogeneously isotropic two-phase model is used to test the effectiveness of the PML absorbing boundary condition. Figure 1a shows the original model. The basic material proper-ties of the model are listed in Table 1. The size of the model is 150×150 nodes with 20 cells of PML on each side of the computational domain. The spatial and temporal steps are Δx = Δz = 10 m and Δt = 0.001 s, respectively. A line source with a dominant frequency of 20 Hz is located at the coordinate (50, 50). The material parameters of the model are listed in Table 1. In order to see clearly the absorbing effect of the PML absorbing boundary condition, we enlarge the original model from 150×150 to 256×256 and retain other pa-rameters unchanged.

Figure  1.  (a) Sketch-map of the original homogeneous model with PML absorbing boundary condition. A and B are the outer and inner boundaries of the PML layer; (b) waveform of vertical component of the solid displacement at one point along a vertical line running through the source. Solid line. result from the model with PML; dashed lines. result with reflection.

DownLoad: Full-Size Img PowerPoint

Table  1.  Material and Biot's parameters for a homogeneous fluid-saturated medium (from Pei, 2006)

 | Show Table

DownLoad: CSV

To examine the reflection from the PML interface, we put the receiver along a vertical line running through the source at the point (50, 161) in the enlarged model. The waveforms of vertical displace-ment in solid are shown in Fig. 1b. In the figure, the reference waveform without reflection that is obtained from a much larger model in which the reflections has not arrived within the time window of interest is also shown to compare with the waveform of the PML.

Figure 2a shows snapshot of the vertical compo-nent of the displacement of the solid at 0.35 s when the dissipation coefficient b has the realistic value of 3.00. Figure 2b shows snapshot of the vertical com-ponent of the displacement of the solid of the enlarged model at 0.40 s.

Figure  2.  Snapshots of the vertical components of the solid displacement at 0.35 s (a) in original model and 0.40 s (b) in enlarged model. The model parameters are listed in Table 1; a line source is used.

DownLoad: Full-Size Img PowerPoint

From the above results, we can see that the PML method provides an effective attenuation to outgoing waves. Numerical experiments also demonstrate that this method is stable so long as when the dissipation coefficient b is small.

A Two-Layer Model

In order to illustrate the behavior of the algorithm proposed in this article for heterogeneous media, we test the PML absorbing boundary condition on a two-layer model. In contrast to the homogeneous model, multiwave phases (such as transmitted waves, reflected waves, and converted waves) will be ob-served. As we all know, some of the above-mentioned wave phases are relatively weak compared to the di-rect wave phases and are vulnerable to the artificial reflections from the boundaries. Therefore, whether these wave arrivals can be distinguished determines the effectiveness of the absorbing boundary condition used. The geometry of the model is shown in Fig. 3. An oblique concentration force source with the domi-nant frequency of 20 Hz is located at the center of the model, which is in the upper layer, and the model properties are shown in Table 2.

Figure  3.  Testing model with a horizontal interface.

DownLoad: Full-Size Img PowerPoint

Table  2.  Material parameters for a two-layer model of fluid-saturated medium (from Dai et al., 1995)

 | Show Table

DownLoad: CSV

The snapshots in Fig. 4a at t=0.20 s shows that the solid and fluid particle displacements are in phase for the fast P-wave and out of phase for slow P-wave, and at the same time, the fast P-wave begins to hit the interface. At time t=0.30 s (in Fig. 4b), transmitted and reflected waves are produced from the fast P-wave, among the transmitted waves are the fast P-waves and the converted S-waves. At time t=0.35 s (in Fig. 4c), the slow P-wave arrives at the interface and produces reflected and transmitted waves. At the same time, the reflected fast P-wave overtakes the slow P-wave and the transmitted fast P-wave begins to hit the boundary of the computational region. Figure 4d shows the wavefield at t=0.40 s, a very small reflection resulting from discretization in the PML scheme, which occurs from the edge.

Figure  4.  Snapshots of the vertical component of the solid displacement at four different times, 0.20 s (a), 0.30 s (b), 0.35 s (c), and 0.40 s (d), from the convolutional differentiator solutions with PML. The model parameters are listed in Table 2. The direct fast a and slow b P-waves can be observed in the solid; the transmitted P- c and P-to-S converted e waves are clearly visible in the solid. The fast d and slow f reflected P-waves are also very clearly visible in the solid. Also, the transmitted slow g P-wave occurred very clearly. On Figs. 4c and 4d, significant refracted waves are also present, as can be remarked as m and n.

DownLoad: Full-Size Img PowerPoint

In order to compare the PML absorbing boundary condition with CE (Clayton and Engquist) absorbing boundary condition, the snapshots from the method with CE absorbing boundary condition and the one from the method with PML absorbing boundary con-dition are shown together in Fig. 5. Clearly, the PML absorbing boundary condition has a far better per-formance.

Figure  5.  Snapshots at t=0.40 s for the two-layer model from the method with PML absorbing condition (a) and for the same model from the method with CE absorbing boundary condition (b). The model parame-ters are listed in Table 2. On the upper and lower sides of the second snapshot, small parasitic reflections (a) and (b) from the sides can be observed.

DownLoad: Full-Size Img PowerPoint

CONCLUSIONS

A convolutional differentiator method combined with the PML absorbing boundary condition is devel-oped for the simulation of seismic wave propagation in fluid-saturated porous media. In the boundary re-gion of the computational domain, the PMLs are used to attenuate the outgoing seismic waves.

From the numerical examples, we can see that the results simulated by the convolutional differenti-ator scheme have little dispersion, which do not inter-fere with the analysis of the normal wavefield. As a result, it will improve the accuracy of velocity analy-sis and migration imaging. The wavefield simulated by the convolutional differentiator scheme is very good. From Fig. 4, we can see that at t=0.20 s the fast P-wave begins to hit the interface, and at t=0.30 s the fast P- wave propagates to the interface and produces transmitted waves and reflections, among the trans-mitted waves are the fast P-waves and the converted S-waves. In Fig. 4c, at time t=0.35 s, the slow P-wave arrives at the interface and produces reflected and transmitted waves. At the same time, the reflected fast P-wave overtakes the slow P-wave and the transmitted fast P-wave begins to hit the boundary of the compu-tational region.

The experiments demonstrate that, on the inter-face in the two-phase media, three kinds of waves can be produced, such as fast and slow P-waves and the S-wave. The fast and slow P-waves can produce the converted waves when meeting with the interface to produce the transmitted and reflected waves. For the reason that the attenuation factor in the second layer is much larger than that in the first layer, the slow P-wave can be absorbed quickly after coming into be-ing, so in the figures the transmitted slow P-wave cannot be seen. The energy of the slow P-wave is lar-ger in the fluid, while the energy of the fast P-wave is relatively obvious in the solid. As we all know, the propagation property of the wave on the interface is the main content in the seismic wave modeling; in ex-ploration seismology, the primary aim is to obtain the reflected waves coming from the interfaces. From Figs. 4d and 5a, we can observe that the boundary reflections are reduced greatly, and the normal reflections are clearly observed, although a very small reflection resulted from discretization in the PML scheme, which occurs from the edge. Furthermore, since the PML boundary condition is directly incorporated into the wave equations, it can be used in other methods.

ACKNOWLEDGMENTS

This study was supported by the National Natural Science Foundation of China (No. 40804008). The author thanks Professor Li Xiaofan (Institute of Geol-ogy and Geophysics, Chinese Academy of Sciences) for his helpful discussion and the anonymous review-ers for the suggestions to improve the manuscript.