This article appeared in a journal published by Elsevier. The attached copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues. Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited. In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier’s archiving and manuscript policies are encouraged to visit: http://www.elsevier.com/authorsrights
Author's personal copy Journal of Applied Geophysics 104 (2014) 75–89
Contents lists available at ScienceDirect
Journal of Applied Geophysics journal homepage: www.elsevier.com/locate/jappgeo
Finite difference modeling of ultrasonic propagation (coda waves) in digital porous cores with un-split convolutional PML and rotated staggered grid Yan Zhang a, Li-Yun Fu a,⁎, Luxin Zhang b,1, Wei Wei a, Xizhu Guan a a b
Key Laboratory of the Earth's Deep Interior, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
a r t i c l e
i n f o
Article history: Received 28 April 2013 Accepted 20 February 2014 Available online 3 March 2014 Keywords: Digital porous core Rotated staggered-grid Numerical modeling Unsplit convolutional perfectly matched layer Ultrasonic coda attenuation
a b s t r a c t Ultrasonic wave propagation in heterogeneous porous cores under laboratory studies is an extremely complex process involved with strong scattering by microscale heterogeneous structures. The resulting coda waves, as an index to measure scattering attenuation, are recorded as continuous waveforms in the tail portion of wavetrains. Because of the contamination of reflections from the side ends and reverberations between the sample surfaces, it is difficult to extract pure coda waves from ultrasonic measurements for the estimation of the P- and S-coda attenuation quality factors. Comparisons of numerical and experimental ultrasonic wave propagation in heterogeneous porous cores can give important insight into understanding the effect of boundary reflections on the P- and S-codas in the laboratory experiment. It challenges numerical modeling techniques by three major issues: the creation of a digital core model to map heterogeneous rock properties in detail, the perfect simulation with a controllable and accurate absorbing boundary, and overcoming the numerical dispersions resulting from high-frequency propagation and strong heterogeneity in material. A rotated staggered-grid finitedifference method of Biot's poroelastic equations is presented with an unsplit convolutional perfectly matched layer (CPML) absorbing boundary to simulate poroelastic wave propagation in isotropic and fluid-saturated porous media. The contamination of boundary reflections on coda waves is controlled by the CPML absorbing coefficients for the comparison between numerical and experimental ultrasonic waveforms. Numerical examples with a digital porous core demonstrate that the boundary reflections contaminate coda waves seriously, causing much larger coda quality factors and thus underestimating scattering attenuation. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Wave scattering at short wavelengths has long been interesting to geophysicists (e.g., Fehler, 1982; Nishizawa et al., 1997; Sato and Fehler, 1998; Wu, 1982). Laboratory ultrasonic measurements are important for such studies where ultrasonic waves interact with small-scale random heterogeneities on a scale of micrometers. However, ultrasonic wave propagation in heterogeneous porous cores is an extremely complex process where scattering effects by individual pores and grains are generally neglected in studying ultrasonic wave propagation. Laboratory measurements have shown that the attenuation level predicted by the combined effect of various absorption mechanisms (e.g., the Biot, mesoscopic-loss, and squirt-flow mechanisms) underestimates the measured level of dispersion and attenuation in rocks (e.g., Arntsen and Carcione, 2001; Dvorkin and Nur, 1995; Mavko et al., 1998). Carcione and Picotti (2006) discussed ⁎ Corresponding author. E-mail addresses:
[email protected] (Y. Zhang),
[email protected] (L.-Y. Fu). 1 Now at Schlumberger Western Geco, Houston, TX, USA.
http://dx.doi.org/10.1016/j.jappgeo.2014.02.012 0926-9851/© 2014 Elsevier B.V. All rights reserved.
the significance of the mesoscopic loss at seismic frequencies and compare it with other mechanisms, such as scattering. Ultrasonic scattering attenuation may occur, particularly when ultrasonic wavelengths are comparable to the scale of pores and grains where the scattering effect will be significant (Wu, 1989). Ultrasonic scattering attenuation by small-scale random heterogeneities in porous cores can be measured by ultrasonic coda waves. However, in laboratory ultrasonic measurements only direct waves are used for velocity and attenuation estimations. The tail portion of an ultrasonic wavetrain is often ignored, because of the sample-size limitation of experiments, the contamination of boundary reflections, the unknown heterogeneity in rocks, and the complexity of received waveforms (Stacey and Gladwin, 1981). The ultrasonic coda as a continuous waveform in the tail portion of an ultrasonic wavetrain is composed of a superposition of incoherent scattered waves by microscale heterogeneities in porous cores. It has been used to measure scattering attenuation in the ultrasonic frequency range (Guo and Fu, 2007; Guo et al., 2009). Because of the contamination of boundary reflections from the side ends of a sample core, however, it is difficult to extract pure coda waves from ultrasonic measurements. Numerical simulations with
Author's personal copy 76
Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
absorbing boundary for ultrasonic wave propagation in digital heterogeneous cores can offer some vital information on how the boundary reflections affect the P- and S-coda waves in laboratory setups. To this end, a rotated staggered-grid finite-difference (FD) method of Biot's poroelastic equations is presented in this article with an unsplit convolutional perfectly matched layer (CPML) absorbing boundary. We conduct comparisons of experimental and numerical wave propagation in heterogeneous porous cores, which can give important insight into understanding the complex process of laboratory ultrasonic propagation. FD numerical implementations of poroelastic wave equations have been extensively carried out over the past decades (e.g., Carcione and Goode, 1995; Carcione and Helle, 1999; Dai et al., 1995; Hassanzadeh, 1991; Masson and Pride, 2007; Masson et al., 2006; Sheen et al., 2006; Wang et al., 2003; Wenzlau and Muller, 2009; Zhu and McMechan, 1991), some of which directly solve Biot's poroelastic equations for the whole model that allow the physical properties to vary both laterally and vertically. A comprehensive review and mathematical details are discussed by Carcione (2007). The conventional PML boundary condition has been successfully implemented in the conventional (Zeng et al., 2001) and staggered-grid (Zhao et al., 2007) FD simulation of poroelastic wave propagation. In general, poroelastic effects are much more pronounced at sonic and ultrasonic frequencies than at seismic frequencies. Gurevich (1996) suggested that all numerical simulations based on complex rheological models should be compared to an equivalent elastic model. This invokes comparisons of poroelastic effects between experimental and numerical data. Chin et al. (1985) used a generalized ray expansion algorithm to analyze the attenuation in Plona's experiments (Plona, 1980). Gurevich et al. (1999) compared experiments on a sample made of sintered glass beads to synthetic seismograms by a global matrix approach. Arntsen and Carcione (2001) simulated the Biot slow wave based on the experimental data (Kelder and Smeulders, 1997) in water-saturated Nivelsteiner Sandstone. More recent poroelastic numerical models focused on the effects of partial saturation (Carcione et al., 2003; Helle et al., 2003; Picotti et al., 2007) and rock heterogeneity (Carcione and Picotti, 2006). Because digital core technology based on X-ray tomography are used increasingly to visualize the fluid distribution and spatial heterogeneities in real rocks, it is possible to simulate poroelastic propagation in authentic heterogeneous rock samples, which may contribute to the interpretation of laboratory test series. There are several major factors challenging numerical modeling techniques for comparisons to laboratory measurements. First, numerical simulations need a digital core model that is expected to reflect most of the characteristics of the true one. Numerical calculations (Arns et al., 2002) on digital micro-tomographic images are the key to map heterogeneous rock properties in detail, including both fluids and grains. Errors in numerical elastic properties of core material will seriously affect comparisons of experimental and numerical propagation in porous cores. Second, to handle the contamination of boundary reflections on the tail portion of ultrasonic wavetrains, numerical schemes need to incorporate a controllable and accurate absorbing boundary algorithm. The classical split PML and conventional unsplit PML methods suffer from large spurious reflections at grazing incidence and low-frequency numerical simulation. An unsplit convolutional PML (CPML) overcomes this difficulty with less memory and more computational efficiency (Komatitsch and Martin, 2007; Martin and Komatitsch, 2009). Martin et al. (2008) incorporated the CPML into the standard staggered-grid FD simulation of Biot's equations to improve absorbing performance at grazing incidence. The last factor, concerning heterogeneities at the pore scale from fluid to solid faces, challenges numerical methods for accurate simulation of subtle transmission/scattering effects across pores and grains in digital cores. Numerical dispersions resulting from high-frequency propagation and strong contrasts in material will easily lead to computational instabilities, if numerical schemes are not designed properly. For instance, the standard staggered-grid FD operators (Virieux, 1986) cause instability problems for strong heterogeneities in media, whereas a rotated staggered grid (RSG) (Saenger and Bohlen, 2004;
Saenger and Shapiro, 2002; Saenger et al., 2000) can improve numerical stabilities for high-contrast discontinuities. Because no averaging of elastic moduli is needed in an elementary cell, the RSG FD scheme can characterize microscales fully by individually discretizing the pore content and the solid matrix (Saenger et al., 2005). Our motivation in this article is to understand the effect of boundary reflections on P- and S-codas in the laboratory experiment based on the comparison between the laboratory records and the numerical waveforms calculated by different thickness of absorbing layers. The CPML and RSG approaches are incorporated into the FD simulation of Biot's equations for ultrasonic wave propagation in a fully saturated digital core to improve stability conditions and absorbing performance at grazing incidence. First, the velocity–stress CPML formulation of Biot's poroelastic equations is presented with CPML parameters setup. Then, we describe the RSG finite-difference implementation with a discussion of stability conditions, and a simulation test is also included for the FD simulation with RSG and CPML. We simplify a 3D heterogeneous core to be a 2D double-phase medium and proceed with poroelastic simulations for ultrasonic propagation. We map a real core for elastic properties based on its digital micro-tomographic image. The transmission of ultrasonic waves through the digital core is simulated with the observation setup exactly the same as laboratory measurements. The numerical waveforms computed by different thickness of absorbing layers are compared with the experimental ultrasonic waveforms. Particular attention is paid to the influence of the boundary reflections on the Pand S-codas. Numerical simulations confirm that the boundary reflections seriously contaminate coda waves by stacking over the tail portion of ultrasonic wavetrains, causing much stronger coda quality factors, and thus underestimating scattering attenuation. 2. Unsplit convolutional PML formulation for velocity–stress poroelastic equations The first-order velocity–stress formulation of Biot's equation is easily reformulated using the classical split PML and conventional unsplit PML methods (described in Appendices A and B). There are two main drawbacks associated with these PML approaches: the requirement of much memory and computation, and more importantly, the poor efficiency at grazing incidence after discretization. The CPML technique (e.g., Komatitsch and Martin, 2007; Roden and Gedney, 2000) has been developed for improvement at grazing incidence and low frequency with less memory and more computational efficiency. This section describes a technique to incorporate the CPML into the velocity–stress formulation of Biot's equation. The technique is developed by introducing memory variables into the CPML model to not have to explicitly store all the past states of the medium during the convolution operation, but rather to calculate this convolution in a recursive way. The key factor to reduce the efficiency of the traditional PML models at grazing incidence is the complex coefficient Sk in Eq. (B3). The CPML approach improves the conventional unsplit PML algorithm by introducing not only the damping profile dk, but also two other real variables ak ≥ 0 and χk ≥ 1, such that Eq. (B3) becomes Sk ¼ χ k þ
dk : ak þ iω
ð1Þ
Substituting Eq. (1) into Eq. (B2) yields (with ∂=∂e xk denoted by ∂e xk ) ∂e xk ¼
" # 1 1 dk 1 ∂xk ¼ − 2 ∂x þ ψk ; ∂xk ¼ Sk χ k χ k ðiω þ ak þ dk =χ k Þ χk k
ð2Þ
where ψk is a memory variable that can be rewritten as d d iωψk þ ak þ k ψk ¼ − k2 ∂xk : χk χk
ð3Þ
Author's personal copy Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
Using an inverse Fourier transform to Eq. (3), we have ∂ψk d d þ ak þ k ψk ¼ − k2 ∂xk ; χk ∂t χk
ð4Þ
where ψk is the inverse of ψk . This equation has an iterative solution of the form: n
n−1 −ðak þdk =χ k ÞΔt
ψk ¼ ψk
e
þ
dk ∂xk −ða þd =χ ÞΔt e k k k −1 ; χ k ðak χ k þ dk Þ
ð5Þ
with the initial condition of iterative solution ψk|t = 0 = 0. Eq. (5) is derived simply by solving a first-order differential equation, leading to the same result in Komatitsch and Martin (2007) obtained by a recursive convolution method. We see that the implementation of the CPML scheme is efficient through introducing the memory variable ψk whose time evolution is governed at each time step by Eq. (5). It avoids splitting the fields in the classical split PML method that involves extensive memory and computation. It also solves the problem of the expensive calculation of convolutions in the conventional unsplit PML approach that requires at each time step a sum over all the previous time steps. The CPML formulation can be implemented in the PML region easily in an existing finite-difference code (without PML) by simply replacing each spatial derivative ∂x with ∂x/χk + ψx and advancing ψk in time using Eq. (5). In this article, the CPML technique based on Eqs. (1) and (5) is used for the first-order formulation of the poroelastic wave equation. Following Collino and Tsogka (2001), the damping profile dk in the PML region is chosen as dk ðr Þ ¼ dmax
r m L
;
ð6Þ
where r is the distance from the calculation point to the inner boundary of the PML region, satisfying 0 ≤ r ≤ L with L being the thickness of the absorbing layer. The constant dmax is set as dmax ¼ −
ðm þ 1ÞV max ln ðRÞ; 2L
ð7Þ
where Vmax is the maximum unrelaxed speed of the pressure wave, m is the order of the multinomial, usually chosen as 2 or 3, and R is the theoretical reflection coefficient, set to 10−6 here. When ak = 0 and χk = 1, Sk reduces to that of the classical PML method. As in Martin and Komatitsch (2009), we make ak and χk vary linearly in the PML layer by taking r ak ðr Þ ¼ πamax 1− ; L
ð8Þ
r m χ k ðr Þ ¼ 1 þ ðχ max −1Þ ; L
ð9Þ
where the maximum value amax = πf0 set at the inner boundary of the PML region with f0 being the dominant frequency of the source wavelet. The maximum value χmax set at the external boundary of the PML region can be determined by a series of simulations.
77
suffer from the problem of numerical dispersion and instability, especially for high-frequency waves propagating in strongly heterogeneous poroelastic media. The rotated staggered-grid differential operators approximate derivatives first along diagonal lines, followed with the derivatives along the Cartesian axes computed by linear mapping. It can improve numerical stabilities for high-contrast discontinuities. In this section, we aim to incorporate the rotated staggered-grid FD algorithm into the CPML implementation for the first-order velocity–stress formulation of poroelastic equations. Fig. 1 compares the stencils of standard and rotated staggered-grid FD operators. The standard method generally stores velocity components, stress components, and physical properties in four grids crosswise. Therefore, it is necessary to interpolate the density and the shear modulus by averaging geometrically over these adjacent grid points for the velocity and stress updates. The averaging is accepted theoretically assuming that (1) the elementary cells are homogeneous, (2) the elastic parameters within an elementary cell are weakly heterogeneous, and (3) the stress components at the surface of the cell are constant. Strong heterogeneities in media will reduce numerical accuracy and cause computational instability easily for high-frequency wave propagation. These problems can be avoided by the rotated approach which, as described in Saenger and Shapiro (2002), defines all the velocity components in one grid and all the stress components and elastic parameters in another grid. For the stress update, the geometrical averaging over these adjacent grid points is not required for shear modulus and other elastic parameters. The arithmetic averaging of densities over the adjacent grid points are only needed for the velocity update. The rotated staggered grid consists of rhombus instead of rectangles, where horizontal or vertical boundaries are jagged. Since the medium parameters are usually discretized in rectangular cells, spatial derivatives need to be rotated for rhombic cells. As introduced by Saenger et al. (2000) for discretization of derivatives in the 2D case, the directions of spatial derivatives for a grid with rectangular elementary cells of length Δx and Δz is rotated from the horizontal/vertical directions x and z to the diagonal directions x and z by 8 Δx Δz > : z ¼ Δx x þ Δz z Δr Δr
ð10Þ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with Δr ¼ Δx2 þ Δz2 . Therefore, the spatial derivatives along the old direction can be rotated in terms of Eq. (10) by a linear combination of the derivatives along the diagonal direction. 8 ∂ Δr ∂ ∂ > > þ < ¼ ∂x 2Δx ∂z ∂x : ∂ Δr ∂ ∂ > > : ¼ − ∂z 2Δz ∂z ∂x
ð11Þ
With Eq. (11), it is easy to define the differentiation operators Dx and Dz that perform the spatial derivatives in x - and z -directions for a frequency-domain field.
3.1. Rotated staggered-grid finite difference method with unsplit convolutional PML
8 1 Δx Δz Δx Δz > > u x− ;z þ ; ω −u x þ ; z− ; ω < Dx uðx; z; ωÞ ¼ Δr 2 2 2 2 : 1 Δx Δz Δx Δz > > : Dz uðx; z; ωÞ ¼ u xþ ;zþ ; ω −u x− ; z− ; ω Δr 2 2 2 2
Standard staggered-grid FD techniques have been widely used for the numerical simulation of wave propagation in general heterogeneous poroelastic media. The standard staggered-grid operators, by approximating derivatives along the horizontal or vertical axis, probably
With Eqs. (11) and (12), the spatial derivatives in x- and z-directions can be rotated to the new grid of rhombic cells by a linear combination of derivatives along the new directions.
3. Rotated staggered-grid FD method with unsplit convolutional PML for poroelastic wave equation
ð12Þ
Author's personal copy 78
Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
Fig. 1. Diagrammatic stencils of standard and rotated staggered-grid FD operators. (a) Standard staggered grid with four grids for velocity, stress and physical properties. (b) Rotated staggered grid with only one grid for all velocity components and one for all stress components and physical properties.
8 > ∂ Δr > < uðx; z; ωÞ≈ ðD uðx; z; ωÞ þ Dz uðx; z; ωÞÞ 2Δz x ∂x : ∂ Δr > > uðx; z; ωÞ≈ : ðDx uðx; z; ωÞ−Dz uðx; z; ωÞÞ 2Δx ∂z
ð13Þ
nλ ¼
A detailed process of the implementation of the rotated staggered grid FD scheme is shown in Appendix C with the CPML absorbing boundary. 3.2. Stability and dispersion Saenger et al. (2000) investigated the numerical stability and grid dispersion of the rotated staggered-grid operator for elastic wave propagation. The resulting stability considerations are valid only in the case of homogeneous media and for second-order operators in time. The stability requirement for elastic-wave propagation holds also for the poroelastic case (Wenzlau and Muller, 2009), but the time step must be small enough for the numerical simulation of the fast wave modes. The stability criterion for the velocity–stress rotated staggered-grid scheme with equal grid spacings in 2D can be expressed as the following inequality n pffiffiffiX ΔtV max ≤C ¼ 1= 2 jck j; Δh k¼1
(minimum wavelength) are required for the rotated staggered-grid eighth-order operator, which is generally calculated by
ð14Þ
where Δt is the time step, Vmax is the maximum wave velocity in the medium, Δh is the spatial step, and ck is the spatial difference coefficients depending on the order of the spatial operator. For an eightorder spatial and second-order temporal operator, C = 0.5497 (Zhang et al., 2010). The accuracy of numerical simulation is governed by numerical dispersion that depends on the number of grid points per smallest wavelength for an order-fixed FD operator. Numerical experiments by Chen et al. (2006) show that no less than 3 grid points per wavelength
V min ; Δhf max
ð15Þ
where Vmin is the minimum wave velocity in the medium and fmax is the maximum frequency of the source (generally fmax = 4f0 with f0 the center frequency of the source). Since dispersion errors accumulate with increasing propagation distance, the sufficiently fine grid spacing Δh is necessary for poroelastic simulation where the slow P-wave is critical for dispersion errors. 3.3. Numerical simulation test We take a double-layer model (800 mm × 400 mm) as example, where the upper layer is an isotropic and homogeneous double-phase medium (470 mm × 400 mm) and the lower is an isotropic elastic homogeneous medium (330 mm × 400 mm). The grid intervals are Δx = Δz = 1.0 mm. Detailed poroelastic parameters are listed in Table 1. The CPML absorbing layer with a thickness of 10 grids is attached around the horizontal and vertical boundaries of the model (see Fig. 2). The solid phase source (Ricker wavelet) is located at x = 400 mm and z = 200 mm, with a center frequency f0 = 120 Hz. The time sampling interval is Δt = 0.0001 ms. The CPML parameters are set as m = 2, R = 10−6, χmax = 1, and amax = 120. These snapshots demonstrate the development of all the phases in the consecutive wavefronts: fast P wave (P1), fast P-wave reflection (P1P1), slow P wave (P2), S-wave converted by fast P wave (P1S), slow P wave converted by fast P wave (P1P2), fast P wave converted by slow P wave (P2P1), and slow P-wave reflection (P2P2). We see a clear onset of all the Pwave circular wavefronts with obvious amplitudes, whereas the converted S wave cannot be observed in the solid-particle velocity
Author's personal copy Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89 Table 1 Parameters of the double layer model. Parameters
Above layer (double phase)
Below layer (single phase)
Κs (GPa) Κd (GPa) μ (GPa) ρs (g/cm3) Ф ν κ (mD) Κf (GPa) ρf (g/cm3) η (MPa)
32 4.8 5.7 2.77 0.3 3 1000 2.25 1 1
12 2.2 2.4 2.25 0.001 2.4 1 – – –
components because of its much smaller amplitude. These wavefront characteristics shown in the snapshots conform to Biot's poroelastic equations and can be compared directly with those from other numerical methods. 4. Laboratory experiment versus numerical simulation 4.1. Rock properties and physical setup of the laboratory ultrasonic experiment The sample used is from the Paleogene formation with depth corresponding to ~2800 m of seismic profiles at a well located off the southern coast of China. The cylindrical core is cut with diameters of 40 mm
79
and lengths of 80 mm. The length-to-diameter ratio of the sample is set to 2 to avoid end constraint effects (Franklin and Dusseault, 1989). Lithologies are homogeneous over the length of the core as determined by visual inspection. The thin section cut from the core is impregnated with colored epoxy to facilitate porosity identification. Petrological examination is carried out on the sample via thin section optical microscopy and mineral identification was augmented by bulk analyses. The sample is a high porosity (19.9%) and moderate permeability (43.2 mD) massive sandstone, with its poroelastic parameters listed in Table 2 for the frame matrix and fluid, respectively. It comprises moderately sorted, subangular to subrounded quartz grains (0.1–0.55 mm in diameter) in point contact with one another. Pore size is variable, ranging from about 0.1 mm up to 0.4 mm. The minor clays that generally reside in pores are mainly comprised of Kaolin clays and fine quartz grains, with the proportion ratio being 2:1 approximately accounting to the image analysis of thin-section microphotographs. Most grain boundaries are direct contacts between rigid framework grains. Intragrain fractures are rare, although grain boundary cracks are ubiquitous. In general, the sample is composed of two parts with different porosities: the solid rock frame (mainly composed of quartz grains) and the clays residing in pores between quartz grains. It should be noted that we use “clays” to indicate the clays residing in pores, which are a loose pack mixed by Kaolin clays and fine quartz grains. Such loosely mixed “clays” generally have higher porosity and permeability than pure clays. The schematic of the physical experimental apparatus is shown in Fig. 3. The acoustic measurements system consists of several key parts. A pulse generator in the PZT-crystal mounted on the steel endplates is
Fig. 2. Wavefield snapshots of all the components for the double-layer model at different times of 0.06 ms and 0.15 ms, respectively. Panels (a) and (b) show solid-particle velocities in the x- and z-directions, respectively; panels (c) and (d) show fluid-particle velocities in the x- and z-directions, respectively.
Author's personal copy 80
Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
Table 2 Parameters of the numerical core model. Parameters
Quartz grain
Clay
Oil
Κs (GPa) Κd (GPa) μs (GPa) ρs (g/cm3) Ф ν κ (mD) Κf (GPa) ρf (g/cm3) η (MPa)
37 24.6 44 2.65 0.05 4.8 10.86 – – –
13.3 1.8 15.6 1.94 0.55 1.16 118.66 – – –
– – – – – – – 0.7 0.88 1
used to generate ultrasonic P and S waves with a 600 kHz center frequency. The receiving transducer is connected to the digitizing board in the PC through a signal amplifier. The amplitudes of transmitted elastic waves are monitored by the digital oscilloscope at the opposite side of the sample. The source pulses last for 0.25 × 10−5 s. The maximum duration of receivers is 0.25 × 10−4 s. The sample is jacketed with a rubber tubing to avoid the shaking between the sample and apparatus. In the laboratory, the rubber jackets are also used to weaken boundary reflections to some degree. The oil-saturated sample is tested under ambient pressure conditions in a triaxial cell along the stress path with a constant effective
pressure of 5 MPa. Under a confining pressure of 10 MPa, the pore pressure increases to 5 MPa by a pore fluid inlet in each endplate allowing passage of pore fluids through the sample. Ultrasonic measurements comprise full waveform records of both transmitted P and S waves at nominal center frequency of 600 kHz. In the laboratory, the receiver and the exciter for P-wave observation are located at the same coordinate on the two sides of the sample. The piezoelectric ceramics in our acoustic receivers are adjusted to record the normal vibrations only when observing P-coda waves, so that most of the transverse motions from P to S waves (i.e. conversed S waves from P waves' scattering) will not be recorded along the normal propagation path. Thus, P-to-S conversion can be neglected as an initial assumption for the estimation of P-coda attenuation by the single-scattering model (Sato, 1977). 4.2. Setup of digital core model How to create a reasonable digital porous core, with a proper assignment of elastic properties for the different components in the model, is a crucial issue for numerical modeling. The rationality of a heterogeneous digital core is always controversial for experimental and numerical comparisons of wave propagation. To make Biot's poroelastic equations applicable for the rotated staggered-grid FD numerical modeling in this study, we need to build a digital core model as a double-phase medium. A 2D digital core model is built in this study based on a digital micrograph of the sample via thin sections that are impregnated with
Fig. 3. Schematic of the experimental apparatus. Labels indicate: (a) transmitting piezoelectric transducer; (b) LVDT: linear variable differential transformer; (c) pore fluid inlet; (d) pore fluid outlet; (e) jacketed rock core sample; (f) radial gauge; (g) ultrasonic platen; (h) load cell; and (i) receiving piezoelectric transducer.
Author's personal copy Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
colored epoxy to facilitate porosity identification. Fig. 4a shows a part of the micrograph at magnification. We see clearly pore geometry and grain size as two major controls on porosity and permeability. It also presents varying amounts of interstitial clays that reduce the averaging pore-throat diameter, therefore affecting permeability. The offwhite color indicates quartz grains and the remainder mainly indicates clays residing in pores. Fig. 4b is the black and white version transformed from Fig. 4a, making it easy to separate quartz grains from clays and also to produce a matrix of two colors (black and white) by RGB identification. The resulting digital core model basically maps heterogeneous rock properties in detail, including both pore and grain structures, in which the white (assigned by 1) indicates quartz grains and the black (assigned by 0) represents interstitial clays residing in pores. Fig. 4c shows the 2D digital core model built for the whole core, with its part magnified into view in Fig. 4d. The whole image, 80 mm × 40 mm in size, has 111,392 pixels (472 × 236) with a resolution of 150 dpi. With such a RGB matrix of the digital core image, the volume ratio of two parts, quartz grains and interstitial clays residing in pores, can be easily calculated as 70%:30% averagely. The RGB identification may not be accurate for the conversion of interstitial clay content, but can mainly capture the heterogeneity scale controlled by pore geometry and grain size which affects scattering significantly.
(a)
81
As described above, the porous core is digitalized in a simple way. Each pixel of the RGB image, however, is a single-phase medium, either the solid (grains/clay) or the fluid (oil). For the implementation of Biot's poroelastic equations using the rotated staggered-grid FD method, each pixel point in the digital core image must become a double-phase medium. Based on the concept of a real reservoir of oil/gas, the fluid exists always as a mixture with sands/clays rather than a pure fluid pool. Therefore, in this study oil is set to be filled in the whole background of the digital core model to assure the model to be a heterogeneous double-phase medium. Each pixel in the digital core image, as a sort of effective double-phase medium, is a mixture of oil with either quartz grains or clays residing in pores. For numerical simulations, we need to define proper physical parameters in each pixel point of the model. Based on the volume permeability (43.2 mD) and porosity (19.9%), the volume ratio (7:3) of quartz grains and interstitial clays (pores), and the proportion ratio (2:1) of Kaolin clays and fine quartz grains in the clays residing in pores, the volume ratio of the grain frame, the total pore space, and the clay content occupying in pores can be calculated as 66.5:20.0:13.5. To simplify the calculation of elastic properties, we use an average modulus estimation (Voigt, 1910), i.e., the arithmetic mean in terms of the proportion of quartz grains and clays. According to Liu and Tang (1998), the porosity of quartz grains generally ranges between 0.1% and 8.7%, by which we
(b)
(d)
(c) Fig. 4. Microphotographs of the sandstone sample tested in the laboratory experiment and the numerical core model built by RGB identification with a resolution of 150 dpi. (a) Close up (at 40 times magnification) of part of the sample. (b) The white and black version of (a) with the white indicating quartz grains and the black mainly representing clays in pores. (c) Numerical core model built with 472 × 236 pixel points for the whole sample (80 mm × 40 mm) where the white and black colors indicating quartz grains and clays, respectively. (d) Close up of part of (c) magnified into view, where the oil saturated originally in pores is assumed to distribute over the whole background so as to make each pixel point as a double-phase medium for numerical simulation. The star and square in (c) denote the positions of the source and receiver, respectively.
Author's personal copy 82
Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
choose a mean value of 5% in this study. Using the volume ratio (7:3) and the measured porosity (19.9%), we can calculate the porosity of the interstitial clays residing in pores as 55%. According to Xue (1986), permeability distribution is mostly related to the porosity of grains and clays. We assume that the ratio of permeability between quartz grains and clays approximately equals to the ratio of their porosities. Using the volume ratio (7:3) and the measured permeability (43.2 mD), we can compute the permeabilities of the quartz grains and the interstitial clays residing in pores as 10.86 mD and 118.66 mD, respectively. As mentioned previously, the interstitial clays residing in pores, are a loose pack mixed by Kaolin clays and fine quartz grains, which generally presents higher porosity and permeability than pure clays. The tortuosity factor, a pure geometry parameter, is computed by ν = 1 − r(1 − 1/ϕ), where r ranges between 0 and 1 for ellipsoid grains (Mavko et al., 1998). With the porosities of quartz grains and clays, we can obtain ν = 4.8 for quartz grains and ν = 1.16 for clays by assuming that r = 1/5. The density of quartz is 2.65 g/cm3 (Mavko et al., 1998). Based on the density of Kaolin clay (1.58 g/cm3) and the volume ratio of Kaolin clays and fine quartz grains for clays, the density of clays occupying in pores can be calculated as 1.94 g/cm3 by the average modulus estimation (Voigt, 1910). Using the bulk moduli 37.0 GPa for quartz grains and 1.5 GPa for Kaolin clays (Mavko et al., 1998), the bulk modulus of clays can be obtained as 13.3 GPa by the average modulus estimation. Similarly, the shear moduli of the core are 44.0 GPa for quartz grains and 15.6 GPa for clays. The dry-rock bulk moduli of the core are computed as 24.6 GPa for quartz grains and 1.8 GPa for clays. These poroelastic parameters are listed in Table 2.
4.3. Setup of computation parameters for numerical simulation To study the influence of boundary reflections on coda waves, we simulate the propagation of poroelastic waves in the heterogeneous digital core model surrounded by CPML layers. We use the ultrasonic source as used in laboratory experiments, which can be represented by a spatial–temporal distribution of body forces (Aki and Richards, 1980) as shown in Fig. 5. The source wavelet has a center frequency of 600 kHz, with the time sampling interval Δt = 0.000008 ms. Because the system of differential equations derives from homogenization considerations, the wavelengths must be larger than a typical averaging elementary volume, which in turn must be at least 10 times larger than the pore or grain sizes (Martin et al., 2008; Pride et al., 2004). This is important numerically to avoid some nonphysical waves and unwanted numerical instabilities or numerical dispersions appearing in the simulations. Based on the ultrasonic measurements of phase velocities, we have Vmax = 5749.9 m/s and Vmin = 3390.6 m/s.
Fig. 5. Source function of pressure used with the time sampling interval Δt = 0.000008 ms and a dominant frequency of 600 kHz.
From the stability conditions (Eqs. (14) and (15)), the spatial and time sampling intervals satisfy Δx = Δz ≤ 0.00047 m and Δt ≤ 0.000045 ms, respectively. In the actual numerical simulation, we take the grid intervals Δx = Δz = 0.00005 m in both the horizontal and vertical directions and the time sampling interval Δt = 0.000008 ms in accordance with that of the source. That is, there are 1600 × 800 grid points for the numerical core model (80 mm × 40 mm). Therefore, we need to interpolate the RGB matrix of the digital core image from 472 × 236 pixels to 1600 × 800 grid points. Numerical simulation is performed with 12,500 time steps for a total duration of 0.0001 s. The point source of pressure is located at the grid point (400, 30) (i.e., x = 20 mm, z = 1.5 mm) in Fig. 4c. In the simulation, the CPML thickness varies from 1 to 30 grids. The CPML parameters are set to be m = 2, R = 10−6, χmax = 1, and amax = 600, 000. To analyze the effect of boundary reflections, the radial and vertical components of the stress and velocity are recorded at the receiver located at the point (20 mm, 75 mm), 100 grid points above the bottom CPML in Fig. 4c. 4.4. Numerical results Fig. 6 shows the resulting wavefield snapshots of solid-particle velocity in the z-directions at different times of 0.0048 ms, 0.0096 ms, 0.0144 ms, 0.0192 ms, 0.024 ms, and 0.0288 ms, with different CPML thicknesses np = 1, 2, 10, 20, and 30 grids, respectively. We see a clear onset with maximum amplitude in consecutive wavefronts as direct ultrasonic waves that will be used to calculate velocity and attenuation. Wave scattering as a superposition of incoherent scattered waves in short wavelengths exhibits strong apparent attenuation due to the small-scale random heterogeneities in the digital core. These snapshots demonstrate the development of ultrasonic scattering waves, as well as the formation of coda as continuous waves in the tail portion of wavetrains. We also see a bundle of repetitive boundary reflections from the surrounding side ends of the numerical core. These boundary reflections, strong at np = 1, become weaker with the increase of CPML thicknesses, almost disappearing at np = 20. The excellent performance of the CPML enables us to control the influence of boundary reflections on coda waves in the digital core model, and further to estimate the contribution of boundary reflections to ultrasonic coda attenuation in laboratory acoustic measurements. Comparisons of numerical and experimental ultrasonic waveforms can provide a major impetus to understand the detailed characteristics of ultrasonic wave propagation in porous core. Particularly, the simulation with absorbing boundary in digital heterogeneous cores offers insight into the effect of boundary reflections on ultrasonic coda attenuation. Fig. 7 compares numerical and experimental records for ultrasonic P- and S-waves. These numerical simulations are conducted with different CPML thicknesses np = 1, 2, 10, 20, and 30 grids, respectively. The coda wave is usually right after the direct wave with its start time marked in the figure. A coda window with the length about 0.02 ms for P coda and 0.017 ms for S coda is selected for both numerical and experimental records. Coda quality factors Qpc and Qsc for both the numerical and experimental ultrasonic records are calculated by the method introduced in Appendix D and marked in the corresponding figures. We see that the boundary reflection is so strong at np = 1 grid (i.e., without the CPML artificial boundary) that significantly contaminates the tail portion of records, and thus makes the codas Qpc and Qsc up to 120.71 and 130.42, respectively, considerably larger than the true values 65.8 and 54.5. In the laboratory, the rubber jackets are usually used around the rock sample in the experiment to weaken boundary reflections to some degree. In addition, the experimental ultrasonic waves in the coda window are low in frequency, as shown in the top panel of Fig. 7, indicating either the scattering attenuation in the high-frequent components or the contamination of boundary reflections by rubber jackets. We note that the numerical simulation slightly lowers the
Author's personal copy Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
83
Fig. 6. Wavefield snapshots of the solid-particle velocity in the z-directions at different times of 0.0048 ms, 0.0096 ms, 0.0144 ms, 0.0192 ms, 0.024 ms, and 0.0288 ms for the numerical core model as shown in Fig. 4c, calculated under different CPML thicknesses np = 1, 2, 10, 20, and 30 grids, respectively.
frequencies of the direct P- and S-waves in contrast to the experimental records, possibly because of the imperfect imaging of poroelastic properties or some theoretical defects in Biot's poroelastic equations. In conclusion, the decay of boundary reflections is fast as increasing the CPML thickness from np = 2 to 10 grids. The resulting codas Qpc and Qsc reduce to a normal level after np = 20 (as shown in Fig. 7) almost without the attribution from boundary reflections. These numerical examples indicate that boundary reflections contaminate coda waves seriously, which can increase the passage of energy in the medium, leading to much larger value of coda quality factor, and decrease the scattering attenuation.
5. Discussions Comparisons of numerical and experimental records for ultrasonic P- and S-waves demonstrate some major departures in waveforms as observed in Fig. 7. The laboratory experiment is carried out on a 3D cylindrical rock sample. Besides the imperfect digital imaging of rock samples, the defective calculation of poroelastic properties in numerical double-phase media, and some theoretical defects in Biot's poroelastic equations, the 3D to 2D simplification of a cylindrical rock sample, assumed in this study, is possibly one of the major reasons for deviations in simulated waveforms. This study focuses on two issues, scattering
Author's personal copy 84
Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
Fig. 7. Comparisons of experimental (top two panels) and numerical (the rest panels) records for ultrasonic P- (left panels) and S-waves (right panels), with the numerical simulation calculated under different CPML thicknesses (marked in the corresponding figure) np = 1, 2, 10, 20, and 30 grids, respectively. The coda window and the calculated coda Qsc value are marked in the figure.
attenuation by small-scale random heterogeneities, and boundary reflections by from the side ends of a sample core. Based on an intuitive understanding, the 3D to 2D simplification may not affect these two issues significantly for the used cylindrical core with its grains and pores being homogeneous in size as determined by visual inspection. The resulting 2D numerical core model mainly captures the scale heterogeneities controlled by pore geometry and grain size which affect scattering significantly. The differences in the boundary reflection between 2D and 3D cases lie mainly in amplitude, which can be determined theoretically. Therefore, comparisons of numerical and experimental records for ultrasonic P- and S-waves, conducted in this study, are reasonable and feasible concerning the two issues to be studied. In this study, Sato's single-scattering model (Sato, 1977) derived from an acoustic wave equation is used for coda attenuation estimations. It assumes that the P-to-S and S-to-P conversions are neglected. Former researches (e.g. Aki, 1992) pointed out that P-to-S scattering is more efficient than S-to-P scattering. As a result, the whole S coda mainly consists of scattered S waves (N90%), whereas S-to-P scattering waves can be neglected (b10%). Therefore, the single-scattering model from acoustic approximation is acceptable for the analysis of S-coda attenuation. However, the same argument may not be true completely for the analysis of Pcoda attenuation because of P-to-S conversion. That is, the whole P coda contains some S waves. To make the single-scattering model applicable
for the estimation of P-coda attenuation, we manage to reduce P-to-S conversion in the laboratory observation. Because the receiver and the exciter are located at the same x-coordinate on the two sides of the sample in the laboratory experiment, the piezoelectric ceramics of the acoustic receiver can be adjusted for P-wave observation to only record normal vibrations, so that most of the transverse vibrations of S waves will not be recorded along the normal propagation path. Thus, P-to-S conversion can be neglected and the whole P coda mainly consists of scattered P waves. In addition, the rock sample used in this study presents a medium heterogeneity, causing weak P-to-S conversion. Knopoff and Hudson (1964, 1967) showed that both P-to-S and S-to-P conversions become weaker for highfrequency propagation. Therefore, the single-scattering model can be able to provide a close estimation for P-coda attenuation (Guo et al., 2009). Generally, the intrinsic quality factor calculated by the spectral ratio method (Toksöz et al., 1979) using direct ultrasonic waves include some scattering loss as the ultrasonic wavelengths (particularly of P waves) are of the same order as the size of grains and pores. Thus, it possibly underestimates the true intrinsic attenuation that depends on the combined effect of various absorption mechanisms (e.g., the Biot, mesoscopic-loss, squirt-flow, and viscous mechanisms). Scattering attenuation depends on small-scale heterogeneities, which, against intrinsic attenuation, may be significant when ultrasonic wavelengths are comparable to the scale of pores and grains. Coda attenuation inferred from the decay rate of
Author's personal copy Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
coda waves (Aki and Chouet, 1975; Singh and Herrmann, 1983) also contains some intrinsic loss. The accurate evaluation of relative contributions of scattering and intrinsic attenuation effects to coda attenuation is still obscure (Herraiz and Espinosa, 1987; Fehler and Sato, 2003). There has been no attempt in this study to evaluate the true intrinsic attenuation using both viscous, Biot, and squirt-flow models, and further to distinguish intrinsic and scattering attenuation mechanisms, which will be studied in the future based on the observed data.
6. Conclusions
Appendix A. Elastic wave equation in saturated poroelastic medium This appendix describes the poroelastic wave equation as the theoretical background for our rotated staggered-grid FD solution with unsplit convolutional PML for simulating ultrasonic wave propagation in digital porous cores. Using tensor notation and neglecting source terms, Biot's equations for an isotropic fluid-saturated porous medium can be written in Cartesian coordinates as (Biot, 1962) 8 ρ u þ ρ f wi ¼ ðλu þ μ Þu j;ji þ μui;jj þ αMw j;ji > < b i ρ f ui þ ρm wi ¼ αMu j;ji þ Mw j;ji −bw˙i ; > : w ¼ ϕ u f −u i
Ultrasonic coda waves, observed as the tail portion of an ultrasonic wavetrain in laboratory ultrasonic measurements, are often ignored as noise because of the contamination of boundary reflections from the side ends of a sample core. Numerical simulations with accuracycontrolled absorbing boundary can provide insight into the effect of boundary reflections on the coda waves in laboratory experiments. Few comprehensive numerical simulations with accurate absorbing boundary have been done for ultrasonic coda waves in digital and heterogeneous porous cores. It really challenges numerical techniques by digital imaging of poroelastic properties, numerical dispersion at high frequency and strong heterogeneity, and accurate absorbing-boundary schemes at grazing incidence. In this article, a rotated staggered-grid finite-difference method of Biot's poroelastic equations is presented to simulate ultrasonic wave propagation in a digital porous core with oil saturated and random heterogeneities. To investigate the influence of boundary reflections on ultrasonic coda waves, an unsplit convolutional perfectly matched layer (CPML) absorbing boundary is employed in the simulation to overcome the difficulty of conventional PML methods at grazing incidence. Numerical experiments with a saturated heterogeneous core demonstrate that the rotated staggered-grid FD scheme with the CPML for ultrasonic wave propagation significantly improve stability conditions at strong heterogeneity and absorbing performance at grazing incidence. The contamination of boundary reflections on coda waves is controlled by the CPML absorbing coefficients for the comparison between numerical and experimental ultrasonic waveforms. We conducted comparisons of numerical and experimental ultrasonic P and S records for a cylindrical rock sample with 40 mm in diameter and 80 mm in length. Numerical simulations with the digital oil-saturated, porous core demonstrate that the boundary reflections contaminate coda waves seriously, causing much stronger coda quality factors and thus underestimating scattering attenuation. The boundary reflections decay fast as increasing CPML thicknesses, almost disappearing at the CPML thickness of 20 grids. The excellent performance of the CPML enables us to estimate the contribution of boundary reflections to ultrasonic coda attenuation in laboratory acoustic measurements. Comparisons of the resulting ultrasonic quality factors of P- and S-codas between the experimental and numerical records indicate that the boundary reflections can increase the passage of energy in the medium, leading to much larger value of coda quality factor, and decrease the scattering attenuation.
85
i
ðA1Þ
i
f
where ui, ui , and wi are, respectively, the solid, fluid, and relative displacements, ϕ is the porosity, the bulk density ρb = ϕρf + (1 − ϕ)ρs, ρf and ρs are, respectively, the fluid and solid densities, the effective fluid density ρm = νρf/ϕ with ν the tortuosity (a parameter characterizing the pore geometry), and the indices i and j can be 1 or 2 here in 2D and with the Einstein convention of implicit summation over a repeated index. As described in Picotti et al. (2007) and Wenzlau and Muller (2009), the following simplified relations for the variables, the poroelastic coefficient of effective stress α, the coupling modulus M, the undrained Lamé parameter λu, and the friction coefficient b of the saturated matrix, can be obtained for an isotropic porous medium α ¼ 1−
Kd ; Ks
ðA2Þ
" #−1 α−ϕ ϕ þ ; Ks Kf
ðA3Þ
2 2 λu ¼ K d − μ þ α M; 3
ðA4Þ
M¼
b¼
1 η ¼ ; k κ
ðA5Þ
with Ks and Kf the bulk moduli of the solid and fluid phases, respectively, Kd and μ are the bulk and shear moduli of the dry porous matrix, respectively, k is the fluid conductivity, η is the fluid viscosity, and κ is the rock permeability. In this article, the classical velocity–stress formulation is employed, an approach proposed for elastic-wave equations by Virieux (1986). The velocity–stress version of Biot's poroelastic Eq. (A1) can be rewritten as a set of four evolution equations of the first order for the field variables, the particle velocities vi ¼ u˙i and qi ¼w˙i, the total stress τij, and the pore pressure p. ρv˙ i ¼ ρm τ ij; j þ ρ f p;i þ ρ f bqi ;
ðA6Þ
Acknowledgments
ρ q˙ i ¼ −ρ f τij; j −ρb p;i −ρb bqi ;
ðA7Þ
We are grateful to Bjarne Almqvist at Uppsala for his constructive review and help in editing the original manuscript. We also thank an anonymous reviewer for his valuable comments and help in editing the original manuscript. We gratefully thank Laboratory of HighTemperature and High-Pressure Geodynamics, Institute of Geology and Geophysics, Chinese Academy of Sciences for great help in preparing experimental data. This research was supported by the National Natural Science Foundation of China (Grant No. 40925013) and the Strategic Leading Science and Technology Programme (Class B) of the Chinese Academy of Sciences (Grant No. XDB10010400).
τ˙ij ¼ μ vi; j þ v j;i þ λu vi;i þ αMqi;i δij ;
ðA8Þ
˙ −M αvi;i þ qi;i ; p¼
ðA9Þ
where ρ = ρmρb − ρfρf. Together with suitable boundary conditions, Eqs. (A6)–(A9) describe poroelastic wave propagation in heterogeneous porous media.
Author's personal copy 86
Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
Several assumptions applied to the Biot linear theory for poroelastic wave propagation (Mavko et al., 1998) may affect numerical simulations of ultrasonic wave propagation in digital porous cores. The small deformation to ensure linear elastic behavior, the statistically isotropic porous medium, the large wavelength compared to microscopic pore scales, and the connected pores considered in the equations can be satisfied for ultrasonic wave propagation in fully saturated high-porosity cores with a small size. It is well known that Biot's poroelastic Eq. (A1) are frequency-dependent, characterized by the Biot critical frequency (typically on the order of megahertz). For the range of high frequencies, the Biot slow P-wave is a propagating wave mode and scattering attenuation may be dominant; whereas in the long-wavelength limit, i.e., at frequencies much lower than the Biot critical frequency, friction between the solid matrix and the pore fluid becomes dominant and this wave mode becomes diffusive. Simulations of the low-frequency Biot poroelastic equations particularly challenge numerical schemes. Appendix B. Conventional PML in velocity–stress poroelastic formulation This appendix describes the conventional PML in the velocity–stress formulation as the base on which the unsplit convolutional PML is developed for simulating ultrasonic wave propagation in digital porous cores. The PML absorbing boundary condition (Berenger, 1994) has proven to be very efficient to minimize spurious reflections at the artificial boundaries of the computational domain. The classic PML has been extensively developed and is now routinely used in both acoustic and elastic problems. It is naturally formulated in terms of the first-order velocity–stress equations in time. To introduce the PML for wave propagation, Eqs. (A6)–(A9) will be reformulated using the complex stretched coordinates within the matched layers surrounding a 2D computational domain. As denoted in Zeng et al. (2001) and Martin et al. (2008), the complex coordinatestretching variables are chosen in the frequency domain as Z i xk e xk ¼ xk − d ðsÞds ðk ¼ 1; 2Þ; ðB1Þ ω 0 k pffiffiffiffiffiffiffiffi where i ¼ −1 and the damping functions dk = 0 inside the computational domain and dk N 0 in the PML region. In the PML formulation, the regular coordinate variables xk is replaced by the complex coordinate variables e xk using the fact that ∂ iω ∂ 1 ∂ ¼ ¼ ; ∂e xk iω þ dk ∂xk Sk ∂xk
ðB2Þ
with Sk ¼ 1 þ
dk : iω
ðB3Þ
The PML formulation is conducted directly by changing original wave equations written in terms of variables xk into new equations written in terms of variables e xk . The classical split PML formulation of the velocity–stress poroelastic Eqs. (A6)–(A9) is implemented by: (i) transforming the system of equations into the frequency-domain form, (ii) mapping the resulting equations from the regular Cartesian coordinates to the complex stretched coordinates, (iii) splitting the velocity and strain fields into two components perpendicular and parallel to absorbing boundaries, and (iv) converting the resulting split equations back to the time domain to obtain the final classical PML formulation of the poroelastic equations in a split form. We take Eq. (A6) in 2D as an example with indices i and j replaced by values 1 and 2 corresponding to coordinates x and z, respectively. The equation is split into (
x
ρv˙x ¼ ρm τxx;x þ ρ f p;x þ ρ f bqx : z ρv˙x ¼ ρm τxz;z
ðB4Þ
By transforming into the frequency domain and reformulating in terms of complex coordinates, Eq. (B4) becomes 8 1 1 x > > < iωρvx ¼ ρm τ xx;x þ ρ f p;x þ ρ f bqx Sx Sx : 1 > z > : iωρvx ¼ ρm τ xz;z Sz
ðB5Þ
Converting back to the time domain by an inverse Fourier transform, we have the classical split PML formulation of Eq. (A6) Z 8 < ρ v˙x þ d vx ¼ ρ τ qx dt x x x m xx;x þ ρ f p;x þ ρ f b qx þ dx : : z z
ρ v˙x þ dz vx ¼ ρm τxz;z
ðB6Þ
This modified wave equation has exponentially decaying plane wave solutions in the PML. As indicated by Komatitsch and Martin (2007), the damping coefficient depends on the direction of wave propagation. It is large for a wave propagating close to normal incidence, but becomes significantly smaller for a wave propagating at grazing incidence. The split fields in the classical PML formulation require extensive memory and computation. The PML model can be derived without splitting the fields, leading to the conventional unsplit PML (Song et al., 2005; Wang and Tang, 2003). Note that the inverse Fourier transform of S−1 is k F
−1
−1 −d t ¼ δðt Þ−dk H ðt Þe k ¼ δðt Þ−ζ k ðt Þ; Sk
ðB7Þ
where δ(t) and H(t) are the Dirac delta and Heaviside functions, respectively. Converting Eq. (B2) back to the time domain, one gets
F
−1
∂ ∂e xk
¼ ð1 þ ζ k ðt ÞÞF
−1
∂ ; ∂xk
ðB8Þ
where * denotes a convolution operation. With Eq. (B8), the unsplit form of Eq. (B6) becomes ρv˙x ¼ ρm ð1 þ ζ x ðt ÞÞτxx;x þ ρm ð1 þ ζ z ðt ÞÞτxz;z þ ρ f ð1 þ ζ x ðt ÞÞp;x þ ρ f bð1 þ ζ x ðt ÞÞqx : ðB9Þ The unsplit formulation simplifies the classical split algorithm without sacrificing the accuracy. It requires nearly the same amount of computer storage as does the split-field approach. Unfortunately, the conventional unsplit PML approach does not give satisfactory results at grazing incidence because the complex coefficient Sk used is the same as the split-field approach. Appendix C. Rotated staggered grid FD scheme combined with CPML absorbing boundary This appendix displays the detailed equations of the implement of rotated staggered grid FD scheme with the CPML absorbing boundary. The rotated staggered-grid FD implementation of the first-order spatial derivatives of the velocity and stress components for solid and liquid phases can be expressed with L-order accuracy as: ( L=2 X ∂vi cn jðx;zÞ ¼ vj 1 1 2Δx i ðxþðn−2Þ;zþðn−2ÞÞ ∂x n¼1
)
þvi jðxþðn−1Þ;z−ðn−1ÞÞ −vi jðx−ðn−1Þ;z−ðn−1ÞÞ −vi jðx−ðn−1Þ;zþðn−1ÞÞ ; 2
2
2
2
2
2
ðC1Þ
Author's personal copy Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
( L=2 X ∂vi cn vj jðx;zÞ ¼ −vi jðxþðn−1Þ;z−ðn−1ÞÞ 1 1 2 ) 2 2Δz i ðxþðn−2Þ;zþðn−2ÞÞ ∂z n¼1
ðC2Þ
2
2
2
(
X c ∂qi n qj j ¼ þ qi jðxþðn−1Þ;z−ðn−1ÞÞ 1 1 2 2 ∂x ðx;zÞ n¼1 2Δx i ðxþðn−2Þ;zþðn−2ÞÞ ) L=2
ðC3Þ
−qi jðx−ðn−1Þ;z−ðn−1ÞÞ −qi jðx−ðn−1Þ;zþðn−1ÞÞ ; 2
2
2
2
( L=2 X ∂qi cn −qi jðxþðn−1Þ;z−ðn−1ÞÞ jðx;zÞ ¼ qj 1 1 2 2 2Δz i ðxþðn−2Þ;zþðn−2ÞÞ ∂z ) n¼1 2
2
1 1 ∂x vx þ ψx;vx þ βM ∂x qx þ ψx;qx ; τ˙xx ¼ ðλu þ2μ Þ χx χx
ðC13Þ
1 1 ∂z vz þ ψz;vz þ βM ∂z qz þ ψz;qz ; τ˙zz ¼ ðλu þ2μ Þ χz χz
ðC14Þ
ðC4Þ
−qi jðx−ðn−1Þ;z−ðn−1ÞÞ þ qi jðx−ðn−1Þ;zþðn−1ÞÞ ; 2
1 1 1 ∂x τ xz þ ψx;τxz þ ∂z τ zz þ ψz;τzz −ρb ∂z p þ ψz;p −ρb bqz ; χx χz χz
ρq˙ z ¼ −ρ f
ðC12Þ
−vi jðx−ðn−1Þ;z−ðn−1ÞÞ þ vi jðx−ðn−1Þ;zþðn−1ÞÞ ; 2
87
τ˙xz ¼ μ
1 1 ∂z vx þ ψz;vx þ ∂x vz þ ψx;vz ; χz χx
ðC15Þ
2
1 1 ∂x vx þ ψx;vx þ ∂z vz þ ψz;vz χx χz 1 1 −M ∂x qx þ ψx;qx þ ∂z qz þ ψz;qz : χx χz
˙ −βM p¼ ∂τij ∂x
jðxþ1;zþ1Þ ¼ 2
2
( L=2 X cn τ j 2Δx ij ðxþn;zþnÞ n¼1
ðC5Þ
)
þτij jðxþn;z−ðn−1ÞÞ −τij jðx−n;z−ðn−1ÞÞ −τ ij jðx−n;zþnÞ ;
∂τij
j 1 1 ¼ ∂z ðxþ2;zþ2Þ
Temporal discretization for Eqs. (C9)–(C16) is based on a secondn + 1/2 for velocity, order staggered scheme with vi((n + 1/2)Δt) = vnþ1=2 i n−1=2 n n i =Δt. p(nΔt) = p for stress and pressure, and Dt ϕ ¼ ϕ −ϕi
L=2 o X cn n τ ij jðxþn;zþnÞ −τ ij jðxþn;z−ðn−1ÞÞ −τ ij jðx−n;z−ðn−1ÞÞ þ τij jðx−n;zþnÞ ; 2Δz n¼1
ðC6Þ
o X c n ∂p n pjðxþn;zþnÞ þ pjðxþn;z−ðn−1ÞÞ −pjðx−n;z−ðn−1ÞÞ −pjðx−n;zþnÞ ; jðxþ1;zþ1Þ ¼ 2 2 2Δx ∂x n¼1 L=2
ðC7Þ
L=2 o X ∂p cn n pjðxþn;zþnÞ −pjðxþn;z−ðn−1ÞÞ −pjðx−n;z−ðn−1ÞÞ þ pjðx−n;zþnÞ ; jðxþ1;zþ1Þ ¼ 2 2 2Δz ∂z n¼1
ðC8Þ where cn is the difference coefficient. A rotated staggered-grid FD scheme with the CPML absorbing boundary for the first-order velocity–stress formulation of 2D poroelastic equations can be implemented with L-order spatial accuracy by (1) applying Eq. (5) to Eqs. (C1)–(C8) for the estimation of the memory variables (ψnx;vi , ψnz;vi , ψnx;qi , ψnz;qi , ψnx;τij , ψnz;τij , ψnx,p, ψnz,p) in the rotated operators; (2) substituting these memory variables into Eq. (2) for the implementation of the CPML model to map the first-order spatial derivatives of all the field components into the complex stretched coordinates; and (3) substituting into the velocity–stress poroelastic Eqs. (A6)–(A9) to yield the final CPML formulation of poroelastic equations. ρv˙x ¼ ρm
1 1 ∂ τ þ ψx;τxx þ ∂ τ þ ψz;τxz χ x x xx χ z z xz
þ ρf
1 ∂ p þ ψx;p χx x
þ ρ f bqx ;
ðC9Þ ρv˙z ¼ ρm
1 1 ∂ τ þ ψx;τxz þ ∂ τ þ ψz;τzz χ x x xz χ z z zz
þ ρf
1 ∂ p þ ψz;p χz z
þ ρ f bqz ;
ðC10Þ ρq˙ x ¼ −ρ f
ðC16Þ
1 1 1 ∂x τxx þ ψx;τxx þ ∂z τxz þ ψz;τxz −ρb ∂x p þ ψx;p −ρb bqx ; χx χz χx
ðC11Þ
Appendix D. Coda measurements by the single-scattering model Coda waves are considered as the result of scattering process caused by small-scale random heterogeneities, providing the decisive information on scattering attenuation. Aki and Chouet (1975) proposed the single-backscattering model for a common source and receiver location to explain the time dependence of scattered energy density. The model has been widely used to estimate the decay rate of coda amplitudes. However, the model assumes that the source of the earthquake and receiver are coincident. As an extension of the single-backscattering model, Sato (1977) formulated the energy density equation by considering the case of non-coincident source and receiver. By assuming the distribution of random scatterers to be homogeneous and isotropic (leading to spherical radiation and isotropic scattering), the energy density for coda waves can be expressed as Eðx; t Þ ¼
W ðωÞg 0 ðωÞ V t −Q −1 ωt K 0 H ðV 0 t−r Þe c ; 2 r 4πr
ðD1Þ
where W is the total energy within the unit frequency band around ω, g0(ω) is the total scattering coefficient, r is the hypocentral distance, V0 is the wave velocity, t is the lapse time measured from the origin, and H is the step function. Function K(ν) is given by Z K ðν Þ ¼
1
dw −1
1 2 1 νþ1 −1 1 ¼ ln ¼ tanh ν ν ν−1 ν 2 −w2 ν
νN1:
ðD2Þ
For the case V0t N r, the coda energy density at frequency f reduces to Eð f jr; t Þ ¼
W ð f Þg 0 ð f Þ −2Q −1 c πft K ðaÞe ; 4πr 2
ðD3Þ
aþ1 and a = t/ts with ts the start time of direct S wave where K ðaÞ ¼ 1a ln a−1 (ts = 0.3 × 10−4 s in this study). Taking into account that energy density is proportional to the mean-square amplitudes of coda waves, we take a natural logarithm of Eq. (D3) and rearrange these terms, leading to ln [A(f|r, t)/K(r, a)] = ln C(f) − (πf/Qc)t, where A(f|r, t) represents the observed root-mean-square (rms) amplitudes of the narrow bandpassfiltered waveforms with a center frequency f, K(r, a) = (1/r)K(a)0.5 and C(f) is a constant. In this study, the frequency range is selected as 525–675 kHz. The attenuation of coda waves can be easily calculated
Author's personal copy 88
Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89
Fig. D-1. Unfiltered P- and S-coda waves (upper panels) of laboratory data and the corresponding functions of ln[A(f|r, t)/k(r, a)] versus lapse time t (lower panels) along with best leastsquare fits of the selected coda windows filtered at center frequency of 600 kHz.
from the slope of a straight line fitting the measured ln[A(f|r, t)/K(r, a)] versus t for a given center frequency (Guo et al., 2009). As an example, Fig. D-1 shows the calculation of coda quality factors Qpc and Qsc from the values of amplitude with time in the coda window by fitting the linear part of ln[A(f|r, t)/K(r, a)]. The direct P and S waves are identified by their travel times. Time lengths for P and S codas are, respectively, about 0.02 ms and 0.017 ms. A band-pass filtering is applied before the calculation of coda quality factors. In order to remove the boundary effects of the data dealing, we use a window of 0.05 ms (3 × 10−5 to 8 × 10−5 s) before Fourier transform, after the band-pass filtering we choose the processed data in accordance with the window of the coda portion to calculate the Qc values (Fig. D-1).
References Aki, K., 1992. Scattering conversions-P to conversion-S versus-S to versus-P. Bull. Seismol. Soc. Am. 82, 1969–1972. Aki, K., Chouet, B., 1975. Origin of coda waves: source, attenuation and scattering effects. J. Geophys. Res. 80, 3322–3342. Aki, K., Richards, P., 1980. Quantitative Seismology, I and II. W. H. Freeman, San Francisco. Arns, C.H., Knackstedt, M.A., Pinczewskiz, W.V., Garboczi, E.J., 2002. Computation of linear elastic properties from microtomographic images: methodology and agreement between theory and experiment. Geophysics 67, 1396–1405. Arntsen, B., Carcione, J.M., 2001. Numerical simulation of the Biot slow wave in watersaturated Nivelsteiner Sandstone. Geophysics 66, 890–896. Berenger, J.P., 1994. A perfectly matched layer for the absorption of electromagnetic waves. J. Comput. Phys. 114 (2), 185–200. Biot, M.A., 1962. Mechanics of deformation and acoustic propagation in porous media. J. Appl. Phys. 33 (4), 1482–1498. Carcione, J.M., 2007. Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media, Second ed. Elsevier Science. Carcione, J.M., Goode, G.Q., 1995. Some aspects of the physics and numerical modeling of Biot compressional waves. J. Comput. Acoust. 3, 261–272. Carcione, J.M., Helle, H.B., 1999. Numerical solution of the poroviscoelastic wave equation on a staggered mesh. J. Comput. Phys. 154, 520–527. Carcione, J.M., Picotti, S., 2006. P-wave seismic attenuation by slow wave diffusion: effects of inhomogeneous rock properties. Geophysics 71, O1–O8. Carcione, J.M., Helle, H.B., Pham, N.H., 2003. White's model for wave propagation in partially saturated rocks: comparison with poroelastic numerical experiments. Geophysics 68, 1389–1398. Chen, H., Wang, X.M., Zhao, H.B., 2006. A rotated staggered grid finite-difference with the absorbing boundary condition of a perfectly matched layer. Chin. Sci. Bull. 51 (17), 1985–1994 (in Chinese). Chin, R.C.Y., Berryman, J.G., Hedstrom, G.W., 1985. Generalized ray expansion for pulse propagation and attenuation in fluid-saturated porous media. Wave Motion 7, 43–65. Collino, F., Tsogka, C., 2001. Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media. Geophysics 66, 294–307. Dai, N., Vafidis, A., Kanasewich, E.R., 1995. Wave propagation in heterogeneous, porous media: a velocity–stress, finite-difference method. Geophysics 60, 327–340. Dvorkin, J.P., Nur, A.M., 1995. Elasticity of high-porosity sandstones: theory for two North Sea datasets. Document ID: 1995–0890. 1995 SEG Annual Meeting, October 8–13, Houston, Texas. Society of Exploration Geophysicists.
Fehler, M., 1982. Interaction of seismic waves with a viscous liquid layer. Bull. Seismol. Soc. Am. 72, 55–72. Fehler, M., Sato, H., 2003. Coda. Pure Appl. Geophys. 160, 541–554. Franklin, J.A., Dusseault, M.B., 1989. Rock Engineering. McGraw-Hill, New York. Guo, M.Q., Fu, L.Y., 2007. Stress associated coda attenuation from ultrasonic waveform measurements. Geophys. Res. Lett. 34, L09307. Guo, M.Q., Fu, L.Y., Ba, J., 2009. Comparison of stress-associated coda attenuation and intrinsic attenuation from ultrasonic measurements. Geophys. J. Int. 178, 447–456. Gurevich, B., 1996. On “Wave propagation in heterogeneous, porous media: a velocity– stress, finite difference method” by Dai, N., Vafidis, A., Kanasewich, E.R. (March– April 1995 Geophysics, 327–340). Geophysics 61, 1230–1232. Gurevich, B., Kelder, O., Smeulders, D.M.J., 1999. Validation of the slow compressional wave in porous media: comparison of experiments and numerical simulations. Transp. Porous Media 36, 149–160. Hassanzadeh, S., 1991. Acoustic modeling in fluid-saturated porous media. Geophysics 56, 424–435. Helle, H.B., Pham, N.H., Carcione, J.M., 2003. Velocity and attenuation in partially saturated rocks: poroelastic numerical experiments. Geophys. Prospect. 51, 551–566. Herraiz, M., Espinosa, A.F., 1987. Coda waves: a review. Pure Appl. Geophys. 125 (4), 499–577. Kelder, O., Smeulders, D.M.J., 1997. Observation of the Biot slow wave in water-saturated Nivelsteiner Sandstone. Geophysics 62, 1794–1796. Knopoff, L., Hudson, J.A., 1964. Scattering of elastic waves by small inhomogeneities. J. Acoust. Soc. Am. 36, 338–343. Knopoff, L., Hudson, J.A., 1967. Frequency dependence of amplitude of scattered elastic waves. J. Acoust. Soc. Am. 42, 18–20. Komatitsch, D., Martin, R., 2007. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation. Geophysics 72, SM155–SM167. Liu, Y.R., Tang, H.M., 1998. Rock Mass Mechanics. Press of China University of Geosciences, Beijing 214 (in Chinese). Martin, R., Komatitsch, D., 2009. An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation. Geophys. J. Int. 179, 333–344. Martin, R., Komatitsch, D., Ezziani, A., 2008. An unsplit convolutional perfectly matched layer improved at grazing incidence for seismic wave equation in poroelastic media. Geophysics 73, T51–T61. Masson, Y.J., Pride, S.R., 2007. Poroelastic finite difference modeling of seismic attenuation and dispersion due to mesoscopic-scale heterogeneity. J. Geophys. Res. 112 (B3), B03204. Masson, Y.J., Pride, S.R., Nihei, K.T., 2006. Finite difference modeling of Biot's poroelastic equations at seismic frequencies. J. Geophys. Res. 111 (B10), 10305. Mavko, G., Mukerji, T., Dvorkin, J., 1998. The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media. Cambridge University Press. Nishizawa, O.T., Lei, S.X., Kuwahara, Y., 1997. Laboratory studies of seismic wave propagation in inhomogeneous media using a laser Doppler vibrometer. Bull. Seismol. Soc. Am. 87, 809–823. Picotti, S., Carcione, J.M., Rubino, J.G., Santos, J.E., 2007. P-wave seismic attenuation by slow-wave diffusion: numerical experiments in partially saturated rocks. Geophysics 72, N11–N21. Plona, T., 1980. Observation of a second bulk compressional wave in a porous medium at ultrasonic frequencies. Appl. Phys. Lett. 36, 259–261. Pride, S.R., Berryman, J.G., Harris, J.M., 2004. Seismic attenuation due to wave-induced flow. J. Geophys. Res. 109, 681–693. Roden, J.A., Gedney, S.D., 2000. Convolution PML (CPML): an efficient FDTD implementation of the CFS-PML for arbitrary media. Microw. Opt. Technol. Lett. 27 (5), 334–339. Saenger, E.H., Bohlen, T., 2004. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid. Geophysics 69, 583–591.
Author's personal copy Y. Zhang et al. / Journal of Applied Geophysics 104 (2014) 75–89 Saenger, E.H., Shapiro, S.A., 2002. Effective velocities in fractured media: a numerical study using the rotated staggered finite difference grid. Geophys. Prospect. 50, 183–194. Saenger, E.H., Gold, N., Shapiro, S.A., 2000. Modeling the propagation of elastic waves using a modified finite-difference grid. Wave Motion 31, 77–92. Saenger, E.H., Shapiro, S.A., Keehm, Y., 2005. Seismic effects of viscous Biot-coupling, finite difference simulations on micro-scale. Geophys. Res. Lett. 32, L14310. Sato, H., 1977. Energy propagation including scattering effect: single isotropic scattering approximation. J. Phys. Earth 25, 27–41. Sato, H., Fehler, M., 1998. Seismic Wave Propagation and Scattering in the Heterogeneous Earth. Springer-Verlag, New York. Sheen, D.H., Tuncay, K., Baag, C.E., Ortoleva, P.J., 2006. Parallel implementation of a velocity–stress staggered-grid finite-differences method for 2D poroelastic wave propagation. Comput. Geosci. 32, 1182–1191. Singh, S.K., Herrmann, R.B., 1983. Regionalization of crustal coda Q in the continental United States. J. Geophys. Res. 88, 527–538. Song, R.L., Ma, J., Wang, K.X., 2005. The application of the nonsplitting perfectly matched layer in numerical modeling of wave propagation in poroelastic media. Appl. Geophys. 2 (4), 216–222. Stacey, G.P., Gladwin, M.T., 1981. Rock mass characterization by velocity and Q measurement with ultrasonics. Anelasticity in the Earth. Geodynamics Series 4. American Geophysical Union, Boulder, CO, pp. 78–82. Sue, Y.Q., 1986. Groundwater Dynamics. Geology Press, Beijing 16 (in Chinese). Toksöz, M.N., Johnson, D.H., Timur, A., 1979. Attenuation of seismic waves in dry and saturated rocks, I: laboratory measurements. Geophysics 44, 681–690.
89
Virieux, J., 1986. P–SV wave propagation in heterogeneous media: velocity stress finitedifference method. Geophysics 51, 889–901. Voigt, W., 1910. Lehrbuch der Kristallphysik. Teubner, Leipzig 954–964 (Reprint 1928). Wang, T., Tang, X.M., 2003. Finite-difference modeling of elastic wave propagation: a nonsplitting perfectly matched layer approach. Geophysics 68 (5), 1749–1755. Wang, X.M., Zhang, H.L., Wang, D., 2003. Modelling of seismic wave propagation in heterogeneous poroelastic media using a high-staggered finite-difference method. Chin. J. Geophys. 46, 842–849 (abstract in English). Wenzlau, F., Muller, T.M., 2009. Finite-difference modeling of wave propagation and diffusion in poroelastic media. Geophysics 74, T55–T66. Wu, R.S., 1982. Attenuation of short period seismic waves due to scattering. Geophys. Res. Lett. 9, 9–12. Wu, R.S., 1989. Seismic wave scattering. In: James, D. (Ed.), The Encyclopedia of Solid Earth Geophysics. Van Nostrand Reinhold, New York, pp. 1166–1187. Zeng, Y.Q., He, J.Q., Liu, Q.H., 2001. The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media. Geophysics 66, 1258–1266. Zhang, L.X., Fu, L.Y., Pei, Z.L., 2010. Finite difference modeling of Biot's poroelastic equations with unsplit convolutional PML and rotated staggered grid. Chi. J. Geophys. (abstract in English) 53, 2470–2483. Zhao, H.B., Wang, X.M., Wang, D., 2007. Application of the boundary absorption using a perfectly matched layer for elastic wave simulation in poroelastic media. Chin. J. Geophys. 50, 581–591 (abstract in English). Zhu, X., McMechan, G.A., 1991. Numerical simulation of seismic responses of poroelastic reservoirs using Biot theory. Geophysics 56, 328–339.
