Poroelastic finite-difference modeling for ultrasonic waves in digital porous cores

1.   Introduction

Coda as the result of scattering processes caused by small-scale heterogeneities has been widely used to measure the small-scale random heterogeneities in the lithosphere (Aki 1969; Aki and Chouet 1975). It is recorded as a continuous wavetrain in the tail portion of seismograms. In the past decades, coda attenuation has been co mmonly measured in the frequency range of 1–30 Hz. It is a useful seismological tool to estimate the strength of heterogeneity in the lithosphere (Sato 1977; Wu and Aki 1985). Wave scattering in short wavelengths has long been interesting to geophysicists. There have been some attempts (e.g., Wu 1982; Fehler 1982) to study scattering attenuation in shorter wavelengths, using seismogram envelopes to reveal information about inhomogeneity on a scale of meters. The envelopes of seismograms of artificial sources by physical modeling have been used to study scattering processes for higher frequencies (Matsunami 1991; Nishizawa et al. 1997; Sivaji et al. 2002; Fukushima et al. 2003). 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. Laboratory measurements have shown that the attenuation level predicted by the combined effect of various mechanisms (e.g., the Biot, mesoscopic-loss, and squirt-flow mechanisms) underestimates the measured level of dispersion and attenuation in rocks (e.g., Dvorkin and Nur 1995; Mavko et al. 1998; Arntsen and Carcione 2001). Carcione and Picotti (2006) discussed the significan ce of the mesoscopic loss at seismic frequencies and compare it with other mechanisms, such as scattering. Ultrasonic scattering attenuation may occur, particularly when wavelengths are comparable to the scale of pores and grains where the scattering effect will be significant (Wu 1989).

Ultrasonic scattering attenuation can be measured by ultrasonic coda waves. The ultrasonic coda as a continuous waveform in the tail portion of ultrasonic wavetrains is composed of a superposition of incoherent scattered waves by microscale heterogeneities in porous cores. However, in many cases, only direct waves in laboratory ultrasonic measurements are used for velocity and attenuation estimations. The tail portion of an ultrasonic wavetrain is often ignored, possibly 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). Guo and Fu (2007) and Guo et al. (2009) made an attempt to measure scattering attenuation in the ultrasonic frequency range using ultrasonic coda waves. 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 absorbing boundary for ultrasonic wave propagation in digital heterogeneous cores can offer crucial information on how the boundary reflections affect the P- and S-codas in laboratory experiments. To this end, a staggered-grid finite-difference (FD) method of Biot's poroelastic equations is presented in this article with unsplit convolutional perfectly matched layer (CPML) absorbing boundary. We conduct comparisons of experimental and numerical wave propagation in heterogeneous porous cores, w hich can give important insight into understanding the complex process of laboratory ultrasonic wave propagation.

Targeted at wave propagation in poroelastic media, various FD numerical modeling techniques for poroelastic wave equations have been extensively studied over the past decades (e.g., Zhu and McMechan 1991; Dai et al. 1995; Carcione and Goode 1995; Carcione and Helle 1999; Wang et al. 2003; Sheen et al. 2006; Masson et al. 2006; Wenzlau and Muller 2009), most of which directly solve Biot's poroelastic equations using conventional staggered-grid FD methods. A comprehensive review and mathematical details can be referred to Carcione (2007). The conventional PML boundary condition has been successfully implemented in the FD (Zeng et al. 2001) and staggered-grid FD (Zhao et al. 2007) 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. Gurevich et al. (1999) compared experiments on a sample made of sintered glass beads to synthetic seismograms by a global matrix approach. Arn tsen 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). Since digital core technology based on X-ray tomography is 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.

Few comprehensive numerical simulations with controllable absorbing boundary have been done for ultrasonic coda waves in digital and heterogeneous porous cores. The simulation of experimental ultrasonic wavetrains allows for the comprehensive investigation of the effect of boundary reflections from the side ends on the P- and S-codas in the laboratory experiment. However, it challenges numerical techniques by three major issues. First, the numerical simulation for comparisons to laboratory measurements needs to abstract a numerical model with its poroelastic properties reflecting the characteristics of true rocks. Micro-tomographic images are made to map heterogeneous rock properties in detail, including both pore and grain structures. Numerical calculations (Arns et al. 2002) on digital image data are the key to obtain accurate data for flow and elastic properties. Improperness of model setting or errors in numerical elastic properties of core material will affect seri ously comparisons of experimental and numerical propagation in porous cores. Second, the strong boundary reflections from the side ends of a sample core contaminate the tail portion of ultrasonic wavetrains seriously, disabling the extraction of coda waves from ultrasonic measurements. 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 CPML overcomes this difficulty with less memory and more computational efficiency (Komatitsch and Martin 2007; Martin and Komatitsch 2009). The CPML has been incorporated into the standard (Martin et al. 2008) and rotated (Zhang et al. 2010) staggered-grid FD simulation of Biot's equations to improve absorbing performance at grazing incidence. The last factor, the strong heterogeneity at the pore scale from flow 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 lead to computational instabilities easily, if numerical schemes are not designed properly. The widely used standard staggered-grid FD operators (Virieux 1986) cause instability problems for strong heterogeneities in medium and high-frequency simulation.

We present a staggered-grid high-order FD method of Biot's equations with its accuracy controllable to simulate the subtle transmission/scattering effects of elastic wave across pores and grains in digital cores. In the first section, we derive the staggered-grid FD formulas with a rbitrary even-order (2L)-order accuracy in space. The high-order FD approximation for spatial derivatives enables a controllable operation to reduce the numerical dispersion for a given accuracy. An alternative way to calculate FD coefficients is provided. The derivation is based on the standard Taylor series expansion but is given in a neat and explicit form. A new layout of grid cells is given to enhance the stability of our formulas. In the second section, we incorporate the high-order staggered FD scheme with the CPML to investigate the effect of boundary reflections on P- and S- codas in digital porous cores. Comparisons between laboratory records and the numerical waveforms with controllable boundary attenuation are made by setting different thicknesses of absorbing layers. In the final section, we use the 2L-order FD scheme with the CPML to synthesize seismograms for digital and heterogeneous porous cores. To avoid the difficulty of 3D digital core image, we simplify the model to be a 2D double-phase medium by digitalizing a section picture of the real core. The transmission of ultrasonic waves through the digital core is simulated with the observation setup exactly the same as laboratory measurements. Numerical experiments with a saturated heterogeneous core demonstrate that the 2L FD scheme with the CPML for ultrasonic wave propagation significantly improves stability conditions at strong heterogeneity and absorbing performance at grazing incidence. Particular attention is paid to the influence of the boundary reflections on the P- and S-codas. Numerical simulations confirm that the boundary reflections contaminate coda waves seriously, which can increase the return of energy into the medium, leading to much larger values of coda quality factor, and decrease the scattering and intrinsic attenuations.

2.   An arbitrary even-order (2L) staggered-grid FD algorithm for velocity-stress equations of porous media

As described in Appendix 1, the first-order velocity-stress formulation of Biot's equation (Eqs. (34)–(37)) is easily solved using the conventional staggered-grid FD techniques. The most attractive feature of the staggered-grid algorithms is that they have higher computational accuracy and better stability than the conventional-grid FD methods. However, the staggered-grid FD algorithms may still suffer from the problem of instability, especially for high-frequency propagation and strong heterogeneity in poroelastic cores. The staggered-grid algorithms have been improved by Mora (1989), Graves (1996), Moczo et al. (2002), Saenger and Bohlen (2004), and others. Based on the standard Taylor series expansion, Pei (2006) derived an explicit staggered-grid FD method for investigation of S-wave splitting and S-wave second splitting at low frequencies in weakly anisotropic media with a low crack density. The method is extended to a high-order 2D 3C FD code for elastic waves propagation in fractured coalbeds with strong anisotropy (Pei et al. 2012). Similarly, Liu and Sen (2009) formulated an implicit staggered-grid high-order FD method by a plane wave theory and the Taylor series expansion. However, t hese works were confined to simulate elastic wave propagation in isotropic and anisotropic media. A staggered-grid high-order FD solution of poroelastic Eqs. (34)–(37), described in this section, can be used to model ultrasonic wave propagation in digital porous cores with strong heterogeneity. Particular attention is paid to the numerical approximation of spatial derivatives.

Assuming L is the half of the differential operator length and u(x) has a (2L + 1)-order derivative, the (2L + 1)-order Taylor expansion of u(x) at with the truncation error O(∆x^2L+3) can be written as

where Δx is the spatial interval in the x-direction and m = 1, 2, …L is the number of staggered-grid points. From Eq. (1), we have

A staggered-grid FD expression of the first-order spatial derivative can be expressed with (2L + 1)-order accuracy as

Substituting Eqs. (2) into (3), we obtain the following 2L-order expression:

where and the truncation error is

Comparing the coefficients of the corresponding terms on the both sides of Eq. (4), we have

The FD coefficients a_m(m = 1, 2, …, L) can be obtained by solving Eq. (6). The solution is

which can be simplified for L→∞ to (Fornberg and Ghrist 1999)

Its coefficients can be also calculated from the staggered convolution differentiator (Zhou et al. 1994) and the Padé expansion (Fornberg and Ghrist 1999), but we derived them in a different way. From Eq. (3), the first-order spatial derivative can be given by

which has 2L-order accuracy.

A staggered-grid FD scheme with 2L-order spatial accuracy and 2nd-order temporal accuracy is implemented for poroelastic wave propagation. All elastic constants and wavefield variables are placed at different positions within a staggered-grid FD cell. The detailed layout of a cell is shown in Fig. 1, and the physical parameters at the corresponding position are listed in Table 1. After a lot of numerical tests, we found that the calculation is a little more stable using the layout shown in Fig. 1 than using other layouts (e.g., Virieux 1986). For the first-order velocity-stress formulation of Biot's equation as described in Appendix 1 (Eqs. 34–37), the discrete form of particle velocity component v_x is given by

Figure  1.  Layout of a finite-difference grid cell in the XZ-plane for the staggered-grid formulation in 2D poroelastic media. The node numbers define the locations of wavefield variables and elastic parameters that are listed in Table 1

DownLoad: Full-Size Img PowerPoint

Table  1.  Distribution of elastic constants and wavefields in a grid cell

 | Show Table

DownLoad: CSV

and

where

f_x is the body force, s = - ϕp, and the discrete forms of stress components are expressed as

and

where ∆t is the time step, and the others are given by

and

with Δz the spatial step in the z-direction.

In the above equations, the spatial-differential functions are approximated by

and

The staggered-grid discrete form of other components (e.g., σ_zz and v_z) can be obtained in the same way.

In numerical examples, we set L = 10 to guarantee a computation accuracy of O(Δt², Δx²⁰). Two steps are needed to perform the scheme for poroelastic wave propagation. First, the particle velocities and stresses are initialized with the initial conditions v_i(x, z, 0) = 0 and σ_ij(x, z, 0) = 0, where the indices i and j represent the spatial coordinate (x, z). Second, the initial model is updated for all time as the following

and

where ∂_j is a short form of the first-order spatial derivative, f_j(x, z, t) the source wavelet in a type of stress, Ω the interior calculation domain, and the absorbing layers to be addressed in the next section. The stability condition of the presented scheme for 2D cases is given by the following inequality:

where V_max is the maximum phase velocity and a_l the FD coefficients.

The accuracy of the proposed high-order staggered FD method is examined by comparing with the exact solution given by Dai et al. (1995) for a homogeneous porous medium with its poroelastic parameters listed in Table 2. The model is 1, 500 × 1, 500 m with the grid intervals being Δx = Δz = 2.5 m. The time sampling interval is Δt = 0.2 ms. The P-wave source function is a Ricker wavelet with a dominant frequency of 50 Hz and is placed in the center of the model. Figure 2 shows an excellent agreement between the numerical (dotted line) and analytical (solid line) solutions at the receiver (x = 750 m, z = 625 m) as functions of time.

Table  2.  Poroelastic parameters for a homogeneous porous medium

 | Show Table

DownLoad: CSV

Figure  2.  Comparison of numerical (dots) and analytical (solid line) solutions of Biot's equation for a homogeneous porous medium with its poroelastic parameters listed in Table 2. The vertical component of the solid phase is shown as function of time. The P-wave source function is a Ricker wavelet with a dominant frequency of 50 Hz

DownLoad: Full-Size Img PowerPoint

3.   Unsplit convolutional PML formulation for velocity-stress poroelastic equations

As described in Appendix 2, the first-order velocity-stress formulation of Biot's equation is easily reformulated using the classical split PML and conventional unsplit PML methods. There are two main drawbacks associated with these PML approaches: the requirement of much memory and computation, and more importantly, the poor performance at grazing incidence after discretization. The CPML technique (e.g., Roden and Gedney 2000; Komatitsch and Martin 2007) has been developed for improvement at grazing incidence and low frequency with less memory and more computational efficiency. This section incorporates the CPML technique 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 leading to the reduction of efficiency of the traditional PML models at grazing incidence is the complex coefficient S_k in Eq. (40). The CPML approach improves the conventional unsplit PML algorithm by introducing not only the damping profile d_k, but also two other real variables a_k≥0 and χ_k≥1, such that Eq. (40) becomes

Substituting Eqs. (18) into (39) yields (with denoted by )

where is a memory variable that can be rewritten as

Using an inverse Fourier transform to Eq. (20), we have

where ψ_k is the inverse of . This equation has an iterative solution of the form:

with the initial condition of iterative solution

Equation (22) 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 recursive Eq. (22). It avoids splitting the fields in the classical split PML method that involves with 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. (22).

In this article, the CPML technique based on Eqs. (18) and (22) is used for the first-order formulation of the poroelastic wave equation. Following Collino and Tsogka (2001), the damping profile d_k in the PML region is chosen as

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 d_max is set as

where V_max 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 a_k and χ_k = 1, S_k reduces to that of the classical PML method. As in Martin and Komatitsch (2009), we make a_k and χ_k vary linearly in the PML layer by taking

where the maximum value a_max = πf₀ set at the inner boundary of the PML region with f₀ 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.

An 2L staggered FD scheme with the CPML absorbing boundary for the first-order velocity-stress formulation of 2D poroelastic equations can be implemented (i) using Eq. (22) for the estimation of the memory variables in the differential operators; (ii) substituting these memory variables into Eq. (19) for the implementation of the CPML model to map the first-order spatial derivatives of all the field components into the complex stretched coordinates; (iii) substituting into the velocity-stress poroelastic Eqs. (34–37) to yield the CPML formulation of poroelastic equations.

A simple two-layer model (314 × 114 m) separated at the depth of 67 m is used to compare the absorbing efficiency between the conventional and unsplit convolutional PML methods, with a particular attention paid to the absorbing performance at grazing incidence. Both the layers are homogeneous porous media with the physical parameters of the model listed in Table 3, showing the same physical parameters in fluid phase, but different in solid phase. The grid and time sampling intervals are Δx = Δz = 0.5 m and Δt = 0.05 ms, respectively. The P-wave source function is a Ricker wavelet with a dominant frequency of 40 Hz and is placed at x = 56 m and z = 97 m. The absorbing layers with a thickness of 10 grids are attached around the vertical and horizontal boundaries of the model. The parameters for the absorbing layers are set to be m = 2 and R = 10⁻⁶. The additional two parameters for the CPML are χ_max = 1.0 and α_max = 40.0. Figure 3 shows snapshots of the simulati on for applying the CPML to all the absorbing boundaries (left panel) and then replacing the bottom boundary by the conventional PML (right panel). We see an excellent performance of the CPML in all the absorbing boundaries, whereas the conventional PML exhibits significant spurious oscillations along the bottom boundary to which the waves reach largely at grazing incidence.

Table  3.  Physical parameters for the two-layer model

 | Show Table

DownLoad: CSV

Figure  3.  Snapshots of the vertical component of the solid phase, at time 0.06 (top), 0.12, 0.18, 0.24, and 0.3 s (bottom), for a two-layer homogeneous porous media with their poroelastic parameters listed in Table 3, where the CPML is implemented on the four sides (left panel), then replaced by the conventional PML only for the bottom boundary (right panel). Amplitudes highe r than a threshold of 5 % of the maximum are shown in red (positive) or blue (negative), with the normalized values raised to the power 0.3 to enhance small amplitudes that otherwise would not be clearly visible. The cross indicates the location of the source, and the interface between the two media is represented by a light line

DownLoad: Full-Size Img PowerPoint

4.   Applications of numerical simulations to digital porous core

The rock sample used in this experiment is a high-porosity (19 %), moderate-permeability (41 mD) massive sandstone. It is from about 3, 300 m depth and presents amounts of intra-grain pores and fractures. Figure 4 gives a close up of the center of the sandstone sample, showing microstructures with various sizes of quartz grains and pores. It comprises moderately sorted, subangular to subrounded quartz grains (0.1–0.3 mm diameter) in point contact with one another. Pore size is variable, ranging from about 0.11 mm up to 0.42 mm. The minor clays and glauconite that are present generally reside in pores, and most grain boundaries are direct contacts between rigid framework grains. Intragrain fractures are rare, although grain boundary cracks are ubiquitous.

Figure  4.  Close up of the center of the sandstone sample tested under triaxial loading in this study, showing microstructures with various sizes of quartz grains and pores

DownLoad: Full-Size Img PowerPoint

How to create a reasonable digital porous core, with a proper assignment of elastic properties for the different components in the model, is always controversial for experimental and numerical comparisons of wave propagation. Numerical calculations (Arns et al. 2002) are conducted to obtain the digital image data for fluid and elastic properties based on the micro-tomographic image of the real core. The resulting digital core model, as shown in Fig. 5, basically maps heterogeneous rock properties in detail, including both pore and grain structures, with the white and black colors indicating the quartz grains and interstitial clays residing in pores, respectively. The crucial problem to be solved is that every point of the image, however, is a single-phase medium, either the quartz grains or the interstitial clays. To make Biot's poroelastic equations applicable for numerical modeling, we need to build a digital core model as a double-phase medium.

Figure  5.  Numerical core model with the white and black colors in dicating the quartz grains and clays, respectively, and the oil saturated in pores is filled in the whole background. The positions of the source and receiver are indicated by star and triangle, respectively. The vertical and horizontal lines around the core represent the CPML layer edge

DownLoad: Full-Size Img PowerPoint

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. That is, every point of the model, as a sort of effective double-phase medium, is a mixture of oil with either quartz grains or interstitial clays residing in pores. It is obvious that the core is significantly simplified with attempt to approximate the main characteristics of the true model. Parameters of the core are shown in Table 4. The rock matrix of the original core consists of two parts: the "harder" rock frame (mainly quartz grains with a little of potash feldspar and glauconite) and the "softer" clays reside in pores between quartz grains (including ankerite, quartz, clay, and opaque). Both have different porosities. The space ratio of the quartz grains (white color) and the clays (black color) residing in pores is about 70 %:30 % based on the image analysis of thin-section microphotographs.

Table  4.  Physical parameters for the numerical core model

 | Show Table

DownLoad: CSV

To compare with the numerical simulation on poroelastic media, we conduct a laboratory experiment using an ultrasonic system. In the Acoustic Measurements System with its schematic diagram shown in Fig. 6, the core is cut to a cylindrical shape, generally 40 mm in diameter and 80 mm in length. It is jacketed with rubber tubing to isolate it from the confining pressure. The length-to-diameter ratio of the rock specimen is set to 2 to avoid end constraint effects (Franklin and Dusseault 1989). The core with oil saturated is tested under ambient pressure condition in a triaxial cell along the stress path with a constant effective pressure. Ultrasonic S waves with a 600-kHz characteristic frequency are emitted from a pulse generator (Panametrics 5077PR) in the PZT-crystal mounted on steel endplates. 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 (Tektronix TDS 420A) at the opposite side of the sample. The source pulses lasted for 0.25 × 10⁻⁵ s. The maximum duration of receiver is 0.25 × 10⁻⁴ s.

Figure  6.  Schematic of the experimental apparatus, consisting of transmitting piezoelectric transducer (indicated by a star), receiving piezoelectric transducer (indicated by a filled triangle), jacketed rock core, rubber jackets, pore fluid inlet, and downstream pore fluid outlet

DownLoad: Full-Size Img PowerPoint

The poroelastic simulation for the heterogeneous core model shown in Fig. 5 is conducted with the accuracy-controllable CPML layers around the model to demonstrate the influence of boundary reflections on coda waves. The ultrasonic S-wave source with a center frequency of 600 kHz used by the laboratory experiment is employed in the numerical example, shifted in time by 0.000008 ms for a null initial condition. The grid intervals are Δx = Δz = 0.00005 m to produce a total size of 1, 600 × 800 grid points (exclude the CPML layers). The time sampling interval is Δt = 0.000008 ms in accordance with the time step of the source function. The thickness of the CPML layers varies from 1 to 30 grid points, with the CPML parameters set to be m = 2, R = 10⁻⁶, χ_max = 1.0, and α_max = 600000.0. The radial and vertical components of the stress and velocity are recorded to analyze the effect of boundary reflections.

Figure 7 shows wavefield snapshots of vertical component at different times for the numerical core model as shown in Fig. 5, with different CPML thicknesses n_p = 1, 5, and 15 grids. We see a clear onset with maximum amplitude in the consecutive wavefronts as the direct ultrasonic waves for the calculation of 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 n_p = 1, become weaker with the increase of CPML thicknesses, almost disappearing at n_p = 15. 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.

Figure  7.  Wavefield snapshots of vertical component at different times of 7.2, 12, 16.8, and 26.4 µs for the numerical core model as shown in Fig. 5, with different CPML thicknesses n_p = 1, 5, and 15 grids

DownLoad: Full-Size Img PowerPoint

Comparisons of numerical and experimental ultrasonic waveforms can provide a major impetus to the understanding of 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. It really challenges the numerical techniques by the digital image of poroelastic properties, the numerical dispersion of waves because of high-frequency propagation and strong heterogeneity, and the accurate absorbing boundary scheme at grazing incidence. Figure 8 compares numerical and laboratory waveforms for an experimental ultrasonic S-wave record. The numerical simulations are conducted with the different CPML thicknesses n_p = 1, 5, and 15 grids. A similar coda window with the time length about 0.017 ms is selected for both numerical and laboratory records. The coda Q_sc values are calculated by Sato's (1977) single scattering model and marked in the figure to demonstrate what a degree the boundary reflections affect ultrasonic coda attenuation.

Figure  8.  Comparisons of laboratory and numerical ultrasonic S-wave trains with the numerical simulation of the core model as shown in Fig. 5 implemented by different CPML thicknesses n_p = 1, 5, and 15 grids, respectively. The coda window and the calculated coda Q_sc value are marked in the figure

DownLoad: Full-Size Img PowerPoint

We see that the boundary reflection is so strong at n_p = 1 grid that significantly contaminates the tail portion of records, and thus makes the coda Q_sc up to 174.3, considerably larger than the laboratory value 52.8. As shown in Fig. 6, the rubber jackets are usually used around the rock sample in the laboratory experiment to weaken boundary reflections to some degree. For more detail, the waves in the coda window become low in frequency, as shown in the top panel of Fig. 8, indicating strong attenuation in high-frequent components by multi-scatterings in the saturated porous media. In the numerical simulation, however, we use the artificial boundary surrounding the digital core as shown in Fig. 6, causing total reflections at np = 0. The faster decay of boundary reflections can be seen at np = 5 and 15 grids in the CPML thickness as shown in the lower two panel of Fig. 8. The resulting coda Q_sc values reduce to 34.2 and 33.3, toward a normal level without the attribution from boundary reflections. This example indicates that the boundary reflection may contribute around one-third of ultrasonic coda attenuation observed in the laboratory experiment.

We note that the FD numerical simulation lowers the frequency of the direct S wave. Theoretically, the high-order FD approximation with much more orders of Taylor series for higher accuracy (e.g., Wang et al. 2003; Liu and Sen 2009; Pei et al. 2012) should not cause such smoothing effect in the resulting synthetic seismograms. It is obviously difficult to investigate the smoothing effect because of the complexity of problem. The main influences on the problem are possibly associated with the imperfect digital imaging of poroelastic properties, the digitalization of porous cores, and the theoretical defect of Biot's poroelastic equations.

5.   Conclusions

Ultrasonic coda waves, observed as the tail portion of an ultrasonic wavetrain in laboratory ultrasonic measurements, are often ignored as noises because of the contamination of boundary reflections from the side ends of a sample core. Numerical simulations with accuracy-controlled 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 ima ge of poroelastic properties, numerical dispersion at high frequency and strong heterogeneity, and accurate absorbing boundary schemes at grazing incidence. In this article, a staggered-grid high-order FD method of Biot's poroelastic equations is presented to model ultrasonic wave propagation in digital porous cores with strong heterogeneity. To investigate the influence of boundary reflections on ultrasonic coda waves, an unsplit convolutional PML 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 2L FD scheme with the CPML for ultrasonic wave propagation significantly improves stability conditions at strong heterogeneity and absorbing performance at grazing incidence. We conduct comparisons of numerical and experimental ultrasonic S waveforms for a cylindrical rock sample with 40 mm in diameter and 80 mm in length. Numerical simulations show that the boundary reflections from the artificial boundary surrounding the digital core decay fast with the increase of CPML thicknesses, almost disappearing at the CPML thickness of 15 grids. The excellent performance of the CPML enables us to estimate the contribution of boundary reflections to ultrasonic coda attenuation in laboratory acoustic measurements. Comparison of the resulting ultrasonic coda Q_sc values between the experimental records and numerical waveforms without the contamination from boundary reflections indicates that the boundary reflection may contribute around one-third of the ultrasonic coda attenuation observed in laboratory experiments.

Acknowledgments We thank anonymous reviewers for their valuab le comments. This research was supported by the National Natural Science Foundation of China (40925013), the Strategic Leading Science and Technology Programme (Class B) of the Chinese Academy of Sciences (Grant No. XDB10010400), and the China National Major Science and Technology Project (2011ZX05023-005-004).

Appendix 1

Elastic wave equation in saturated poroelastic medium

This appendix describes the poroelastic wave equation as the theoretical background for our 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)

where u_i, u_i^f, and w_i are, respectively, the solid, fluid, and relative displacements, is the porosity, the bulk density ρ_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 Biot coefficient of poroelastic medium α, 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

with K_s and K_f the bulk moduli of the solid and fluid phases, respectively, K_d and μ the bulk and shear moduli of the dry porous matrix, respectively, k the fluid conductivity, η the fluid viscosity, and κ 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 Eqs. (27) can be rewritten as a set of four evolution equations of the first order for the field variables, the particle velocities and , the total stress τ_ij, and the pore pressure p.

where ρ = ρ_mρ_b - ρ_fρ_f. Together with suitable boundary conditions, Eqs. (34)–(37) describe poroelastic wave propagation in heterogeneous porous media.

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. (27) 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 2

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 di gital porous cores. The PML absorbing boundary condition (Bérenger 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. (34)–(37) 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 coordinate-stretching variables are chosen in the frequency domain as

where and the damping functions d_k = 0 inside the computational domain and d_k > 0 in the PML region. In the PML formulation, the regular coordinate variables x_k are replaced by the complex coordinate variables using the fact that

with

The PML formulation is conducted directly by changing original wave equations written in terms of variables x_k into new equations written in terms of variables .

The classical split PML formulation of the velocity-stress poroelastic Eqs. (34)–(37) 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. (34) 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

By transforming into the frequency domain and reformulating in terms of complex coordinates, Eq. (41) becomes

Converting back to the time domain by an inverse Fourier transform, we have the classical split PML formulation of Eq. (34)

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 (Wang and Tang 2003; Song et al. 2005). Note that the inverse Fourier transform of S_k⁻¹ is

where δ(t) and H(t) are the Dirac delta and Heaviside functions, respectively. Converting Eq. (39) back to the time domain, one gets

where * denotes a convolution operation. With Eq. (45), the unsplit form of Eq. (43) becomes

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 S_k used is the same as the split-field approach.