Ultrasonic velocity and attenuation during CO2 injection into water-saturated porous sandstone: Measurements using difference seismic tomography

Physics of the Earth and Planetary Interiors 176 (2009) 224–234

Contents lists available at ScienceDirect

Physics of the Earth and Planetary Interiors journal homepage: www.elsevier.com/locate/pepi

Ultrasonic velocity and attenuation during CO2 injection into water-saturated porous sandstone: Measurements using difference seismic tomography Xinglin Lei a,∗ , Ziqiu Xue b,1 a b

Geological Survey of Japan, National Institute of Advanced Industrial Science and Technology (AIST), Central 7, Higashi 1-1, Tsukuba City, Ibaraki 305-8567, Japan Research Institute of Innovative Technology for the Earth, 9-2 Kizugawadai, Kizu-cho, Soraku-gun, Kyoto 619-0292, Japan

a r t i c l e

i n f o

Article history: Received 14 August 2007 Received in revised form 9 April 2009 Accepted 21 June 2009 Keywords: CO2 storage Difference seismic tomography Velocity Attenuation Ultrasonic frequency Laboratory

a b s t r a c t We undertook laboratory-based seismic measurements with dense sensor array at ultrasonic frequencies during the injection of CO2 into a water-saturated sandstone specimen. The resulting high-quality seismic data enabled detailed determination of the relative velocity and attenuation coefficient of the compressional wave using difference seismic tomography, which directly inverses time-lapse changes in rock properties from time-lapse changes in observed data. CO2 migration and water displacement were clearly mapped using tomographic images of relative velocity and the attenuation coefficient. The final and largely stabilised volume fraction of CO2 in the pore space of the sample is about 30–40%. On average, the P-velocity fell by 7.5, 12, and 14.5% and the attenuation coefficient Q−1 increased by factors of 3.3, 2.7, and 3.7 as a result of the replacement of water with CO2 during the injection of gaseous, liquid, and supercritical CO2 , respectively. As a function of gas saturation, both the velocity and attenuation data are in good agreement with results obtained using the White and Dutta–Odé model for partial saturation, indicating that viscous losses due to fluid diffusion are of significant importance for compressional waves travelling at ultrasonic frequencies in porous rocks. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Global warming arising from increasing emissions of greenhouse gases has generally been accepted as a real environmental problem (IPCC, 2001). Carbon dioxide (CO2 ) is the main compound identified as affecting the stability of the Earth’s climate (IPCC, 1995, 2001). To simultaneously stabilize the Earth’s climate and enable the continued use of fossil fuels, thereby avoiding serious harm to the global economy, CO2 capture and storage (CCS) has been proposed as a key technique. The major problem in CCS is storing captured CO2 away from the atmosphere for a sufficiently long time; that is, in the order of 1000 years. The sequestration of CO2 into geological formations has been considered as a good option to solve this problem (Benson, 2005); however, various requirements must be met before the technology will be accepted by the public for wide-scale implementation, including predictions of storage performance, monitoring, verification, and the environmental safety of CO2 storage facilities. Our understanding of all of these issues is dependent upon estimations of the behaviour of CO2 after injection into geological formations, as injected CO2 (whether gas, liquid, or supercritical), does not remain stationary in the sub-

∗ Corresponding author. Tel.: +81 29 861 2468; fax: +81 29 861 3697. E-mail addresses: [email protected] (X. Lei), [email protected] (Z. Xue). 1 Tel.: +81 774 75 2312; fax: +81 774 75 2313. 0031-9201/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.pepi.2009.06.001

surface: it migrates at speeds related to the physical and chemical properties of the CO2 –water system and the geological structures within the formations. Accurate monitoring is necessary to confirm the containment of CO2 , to assess leakage paths, and to improve understanding on interactions between CO2 , pore fluids, and the rock-forming minerals (Benson, 2005). Monitoring is also necessary to quantify the net quantity of CO2 that has been sequestrated within the reservoirs (Chadwick et al., 2004). Seismic surveys currently provide the most attractive approach in obtaining the spatial coverage required for mapping the location and movement of subsurface CO2 (Arts et al., 2004; Li, 2003; Saito et al., 2006; Xue and Lei, 2006; Harris et al., 1995; Chadwick et al., 2004). In general, rock physics applied to fluid substitution problems deals mainly with P-velocity, S-velocity, and density, which are the parameters that are possible to obtain from field data. The seismic wave velocity is a sensitive indicator of lithology, porosity, and saturation. Other rock properties such as the microstructure of pores, permeability, and fluid viscosity play an important role in governing the replacement of a pair of immiscible fluids such as CO2 and water. The replacement of water by CO2 will not only have a dramatic effect on the seismic velocity but also act to reduce the amplitude of seismic waves via anelastic dissipation. Several damping mechanisms have been proposed, including the scattering of small heterogeneities such as cracks, and inertial friction due to viscous flow. In the case of CO2 sequestration, viscous fluid mechanisms might well be the most important factor.

X. Lei, Z. Xue / Physics of the Earth and Planetary Interiors 176 (2009) 224–234

Several approaches based on idealized assumptions (White, 1975; Walsh, 1995) or complex computing (e.g., Tuncay and Corapcioglu, 1996; Lo and Sposito, 2005; Lu and Hanyga, 2005) have been proposed to deal with the inertial friction in porous rocks saturated by two or more immiscible fluids. Most of these works are limited to special models, or the derived equations are too complicated to clarify their physical significance. White (1975) proposed an approximate theory for calculation of the attenuation and dispersion of compressional waves in porous rocks filled mostly with water but containing gas-filled pockets. Several years later, Dutta and Odé (1979a, 1979b) presented a more rigorous solution for White’s model using Biot’s (1962) systematic equations for poroelasticity. The corrected White model is sometimes referred to as the White and Dutta–Odé (hereafter, WDO) model. In the WDO model, velocity and attenuation are described as a function of wave frequency, gas saturation, and the distance between gas pockets in an otherwise water-filled rock. High attenuation at low gas saturation is predicted in the WDO model (Dutta et al., 1979); this has been observed at sonic frequencies in Massilon sandstone (Murphy, 1982). Based on this existing theory, Carcione et al. (2006) presented a methodology used to compute synthetic seismograms for reservoirs subjected to CO2 sequestration. Modelling the seismic properties of reservoir rocks saturated with CO2 and water is a key problem that must be addressed using experimental data; however, little experimental work has been performed in this regard, and seismic attenuation in typical porous rocks containing CO2 and water is not well understood and has yet to be accurately measured. We began an experimental study with the aim of understanding the physical mechanisms governing CO2 –water replacement and improving the accuracy of seismic techniques in monitoring

225

CO2 migration after injection into geological formations. We carried out systematic CO2 -injection experiments in the laboratory using porous rock samples containing inhomogeneous structures and carried out seismic measurements at ultrasonic frequencies with dense sensor array. The obtained high-quality seismic data enabled detailed determination of the relative velocity and attenuation coefficient using difference seismic tomography. We presented preliminary results related to P-velocity in a previous paper (Xue and Lei, 2006), and in the present study will expand upon these results. 2. Experiment and data processing procedure 2.1. Description of experimental setup, test sample and seismic measurement Fig. 1(a) is a simplified diagram of the experiment setup used in this study. The CO2 was injected from the bottom end of the sample at a constant pressure controlled by syringe pump A. Pore pressure was monitored at both ends of the sample. The pore pressure (back pressure) at the top end was kept constant by syringe pump B. Syringe pump C was used to keep the hydrostatic (confining) pressure constant in the vessel. By keeping constant pore pressure at both ends, and holding constant hydrostatic pressure, velocity changes caused by pore pressure build-up were minimised. Thus changes in observed data reflect only displacement between CO2 and water within the test sample. Test samples were shaped into a right circular cylinder with a length range between 100 and 150 mm and a diameter of 50 mm. Stainless-steel end pieces were attached to both sample ends and a grooved plate was placed between the end piece and the sample

Fig. 1. (a) A simplified schematic diagram of the experimental setup. Photographs (b, c) were scanned from a cut surface after the experiment. In (c), pore space is filled with blue epoxy, A and B mark regions of large porosity while C marks a region of low porosity. (d) Ray paths for velocity measurements and cell geometry for tomography. Source and receiver positions are indicated by S# and R#, respectively. (e) A checkerboard model of 4 × 2 cells. (f) Recovered checkerboard. Thick lines show the well resolved area.

226

X. Lei, Z. Xue / Physics of the Earth and Planetary Interiors 176 (2009) 224–234

Fig. 2. Examples of waveform and power spectrum density obtained under dry, wet (fully water-saturated), and partially saturated conditions. The dominant frequency under dry and wet conditions is ∼330 kHz, while that under partial saturation is around 240 kHz.

at each end. Concentric grooves connected by radial grooves were carved on the plate surface to ensure uniform injection over the entire surface of the sample end. The sample assembly was sealed with silicone sealant to prevent immersion of oil, which was used as the hydrostatic pressure medium. We have carried out a systematic experiment using rock cores drilled from Tako sandstone, which was chosen for its high porosity (∼24% under atmosphere pressure). Tako sandstone is finegrained and well-cemented. Measurement using mercury-injection method shows that this rock has a pore-size distribution in the range of 0.01–20 ␮m, with a sharp peak at ∼2 ␮m (Lin et al., 1997; Xue and Ohsumi, 2004a). In the present study, we focused on results obtained from a sample drilled along the bedding plane of the rock. In such a sample, the predominant flow direction of fluids during CO2 injection was along the bedding plane, which possesses a significantly higher permeability than that in other directions. This is generally the case in most field situations. The sample was drilled and cut into a right circular cylinder with a length of 130 mm and a diameter of 50 mm. To directly observe the porosity distribution and other features within the sample, the sample was cut along the axis after the experiment. As shown in Fig. 1(b) and (c), the sample is characterized by strong bedding structures and a non-uniform distribution of porosity. High-porosity regions are clearly identified from the photograph of the surface filled with blue epoxy (Fig. 1(c)). In utilizing the difference tomography introduced above, 16 piezo-electrical transducers (PZTs) were mounted on the sample surface in two parallel lines along the axial direction. The sample was jacked using a carefully painted silicon rubber layer of ∼10 mm thickness and then set in a pressure vessel. Hydraulic oil was used as the confining medium. Fig. 1(d) shows the positions of sources (S1–S8) and receivers (R1–R8), ray paths obtained from ray-tracing using a homogeneous velocity model, and cell geometry used for tomography. The sources were successively connected to a pulse generator through a fast-switching system. The signal received at each receiver was fed to A/D converting equipment and digitized with a dynamic range of 14 bits and a sampling rate of 200 MHz. To enhance the ratio of signal to noise, the waveform for every path was stacked by 64 times. Seismic measurements were performed at a time interval of ∼2 min during CO2 injections. This is the minimum interval of the system used in the present work. We have since developed a new system and shorten the minimum interval to several seconds.

An ultrasonic signal was transformed through 16 feed-throughs that were shielded and absolute with each other in order to minimize background noise and cross-talk. The PZT sensor has a natural frequency of 1 MHz, and the resulting dominant frequency is about 250–350 kHz (Fig. 2), which is far from the natural frequency. In addition, six cross-type strain gauges were mounted on the surface of the sample to measure local strains. 2.2. Experimental procedure 2.2.1. Dry compression The bulk and shear moduli of a dry sample are basic properties in poroelasticity. Compression effect on porosity is also an important factor for data interpretation. Hence, we first carried out a hydrostatic compression test. The sample was dried at room conditions for several weeks and further vacuumed for several hours after being placed in the vessel. Once the dry state was verified, the confining pressure was slowly increased to ∼10 MPa. 2.2.2. Water saturation Next, fresh distilled water was injected from the top end of the sample at a fluid pressure of 1 MPa while the bottom side was closed. This procedure takes more than 2 days to ensure full saturation. 2.2.3. CO2 injection Gaseous CO2 was firstly injected from the bottom side at a pressure of 5 MPa while the pressure at the outlet (upper side) was kept constant at 3 MPa. Three successive runs of liquid CO2 were injected with inlet–outlet pressures of 8–6 MPa, 10–8 MPa, and 12–10 MPa. Finally, supercritical CO2 was injected at inlet–outlet pressures of 12–10 MPa. Before every successive injection, imbibition with fresh distilled water was continuously performed for more than 2 days to ensure that the sample was fully saturated with water. During every run, seismic measurement was operated with a time interval of ∼2 min. Injection pressure, pore pressure at the outlet side, confining pressure, flow rate at inlet and outlet, fluid volume in all pump tanks were digitized and recorded. Experiments were performed in an air-conditioned room and the room temperature was maintained at 23 ◦ C. However, during the injection of supercritical CO2 , the pressure vessel was heated to 34 ◦ C by a controlled heating system. The syringe pump tanks

X. Lei, Z. Xue / Physics of the Earth and Planetary Interiors 176 (2009) 224–234

227

Table 1 Final changes in Vp and attenuation with the injection of different states of CO2 .

Gaseous CO2 Liquid CO2 Supercritical CO2

Pin (MPa)

Pout (MPa)

Pc (MPa)

dVp/Vp(0) (%)

dA (1/m)

dQ−1 /Q−1 (0) (%)

5 12 12

3 10 10

10 15 15

7.5 12 14.5

20 18 23

330 270 370

A, B, and C, containing CO2 , water, and confining oil, respectively, were warmed by circulating water at 34 ◦ C. Confining pressure was controlled independently of pore pressure and was kept constant at 10 MPa or 15 MPa during the experiments. The experimental conditions are summarized in Table 1. Under these conditions, the phase of injected CO2 would keep no change. The physical properties of water and CO2 would be approximately constant during the injection tests. In the present paper, we focused on results associated with the injection of CO2 . Results concerning imbibition are also of interest; these will be presented in a later paper.

ε=

Vmax − Vmin , Vmin

(5)

where  is the angle between the wave vector (ray path) and the axis of symmetry (direction normal to the bedding plane), ε express the degree of anisotropy, Vmax = V(90), and Vmin = V(0). The constant ε describes the fractional difference of the P-wave velocity in the vertical (normal to the bedding plane) and horizontal (parallel with the bedding plane) directions and therefore best describes what is usually called the P-wave anisotropy. Indeed (5) is a special case (ı = ε) of the Thomsen’s notation: V () = Vmin (1 + ı sin2  cos2  + ε sin4 )

2.3. Difference seismic tomography In order to precisely determine relative changes in velocity and attenuation coefficient caused by CO2 replacing water, we developed a two-dimensional difference tomography method. The method used in difference velocity tomography is presented in detail in a previous paper (Xue and Lei, 2006); however, we briefly describe the method here for the sake of completeness. The travel time Tij of a body-wave along the ray path (lij ) from the ith source to the jth receiver is expressed as a path integral of the slowness field p (=1/c, where c is the phase velocity):



Tij =

p ds.

(1)

lij

The travel time difference Tij due to slowness change p is written as



Tij = Tij − Tij =



p ds − lij

l

 ij

(p + p) ds,

(2)

where lij is the new ray path corresponding to the disturbed slowness field. If p is small, change in the ray path due to p can be ignored; thus, Tij can be expressed approximately by the following formula:



Tij ≈

p ds

(3)

lij

Using (3), p can be determined inversely from the differential travel times T with a simultaneous iterative reconstruction technique (SIRT) similar to the classic travel time tomography method based on (1). It is clear that the effect of errors in the slowness p will be partly cancelled out. In addition, differential travel times can be obtained with improved precision via the cross-correlation of waveforms. Therefore, the reconstructed p distribution is expected to be much more accurate than the absolute slowness reconstructed on the base of (1). Moreover, for monitoring fluid migration, we are interested in p much more than p. In general, porous sediments show anisotropy in mechanical properties due to the presence of bedding structures. Thus, velocity and attenuation may depend on the direction of wave propagation. Given that it is reasonable to assume isotropy along the bedding plane, a transverse isotropic model was used in the present study. Further, we used the simplified approximation of P-velocity given by V () = Vmin (1 + ε sin2 ),

(4)

(6)

for weak elastic anisotropy (Thomsen, 1986). However, for our imaging purpose, (4) shows good accuracy. In general, the degree of anisotropy is a function of saturation and may change spatially, thus it is necessary to determine both Vmin (or Vmax ) and ε. However, in such case, the number of unknown parameters would be doubled. In order to reconstruct reliable velocity field with limited observation conditions, we assume a homogenous and constant ε to simplify the inversion problem. The optimum value of ε was estimated for each injection according to the minimum of the standard residual errors of the tomography (Xue and Lei, 2006). Given that we normally inverse Vmin , every result in the present paper relates to Vmin unless stated otherwise. From (4), if we ignore the change in ε we have V ()/V () = Vmin /Vmin , meaning that relative velocity is independent of the direction of wave propagation. Having described the method for difference velocity tomography, we now turn to the method for difference attenuation tomography. The amplitude A of an elastic wave, having a (angular) frequency of ω and phase velocity c, can be expressed as an exponentially decaying function of propagated distance x and the dimensionless quantity Q (Aki and Richards, 2002): A(x) = A0 e−ωx/2cQ ,

(7)

where A0 corresponds to the amplitude of the source at x = 0. In logarithmic form, (7) becomes ln

A(x) = −˛x A0



˛=



ω ω = pQ −1 , 2cQ 2

(8)

where ˛ is the attenuation coefficient. Q−1 is also referred to as the attenuation coefficient in some studies. To avoid confusion, in the present paper we term ˛ and Q the attenuation coefficient and damping factor, respectively. For spatially varying ˛, the body-wave amplitude Aij along the ray path (lij ) from the ith source to the jth receiver can then be expressed as a path integral: ln

Aij A0



=−

˛ ds.

(9)

lij

The relative amplitude, defined as aij = ln(A/Aref ) due to ˛ (change in attenuation coefficient), can thus be written as



aij = aij − aij = −

˛ ds.

(10)

lij

Based on (10), we can determine ˛ inversely from the relative amplitude data in the same way as that used in determining relative velocity.

228

X. Lei, Z. Xue / Physics of the Earth and Planetary Interiors 176 (2009) 224–234

The absolute attenuation coefficient or Q cannot be directly deduced because of the unknown nature of the source amplitude. In fact, the absolute Q is not required if we are only interested in the change in attenuation that arises from the replacement of two fluids; however, if the distribution of Q or its average value can be obtained for either a dry or wet sample, we are then able to convert the relative attenuation coefficient to absolute Q. We need a mean value of absolute Q in comparing our experimental data with theoretical results. For this purpose, we used an aluminium sample having the same geometry and sensor distribution as the test sample to perform reference measurements. We can then use the obtained waveforms as a reference in calibrating the tomographic system and the distribution of Q within the rock sample. Based on Eq. (8), for a given ray path the ratio of the Fourier amplitudes of the sample (A) and the reference (AR ) is ln

x A = (pR QR −1 − pQ −1 )ω AR 2

(11)

Hence, a mean value of Q along the ray path can be determined using the spectral ratio method. Since x is the ray length, (pR QR −1 − pQ −1 ) can be found from the slope of the line fitted to the ratio of amplitude spectrum versus frequency in logarithm coordinates. The damping factor of aluminium QR is very high (>105 ), pQ−1 and thus Q−1 can be determined directly from the slope. In general, we determine a mean value of Q for the water-saturated sample, and then estimate the change in Q–1 from reconstructed average ˛ and p or c using one of the following formulas: Q −1 Q −1

2(p˛ − ˛p) = ωp(p + p) = 2(c˛ + ˛c + c˛)/ω

as the rock sample as described in the previous section. As a result, an average value of Q = 10 was obtained for the water-saturated sample. Second, the differential travel times and relative amplitudes of each measurement were determined from the shift and slope of a linear least-square fit to a target waveform versus the reference waveform. Finally, the time-lapse change in velocity and attenuation was determined relative to the reference data (respect to watersaturated condition) using the difference tomography method described above. Every measurement can be selected as a new reference for subsequent measurements; in fact, difference tomography can be performed for any given pair of measurements. For example, the difference tomography of every measurement made during the injection of CO2 using the standard reference under water-saturated conditions can be used to image the cumulative changes in velocity and attenuation. The relative data between two successive measurements is helpful in understanding the detailed CO2 behaviour in the water-saturated sample. Because the confining pressure and fluid pressures at the inlet and outlet were kept constant during each injection, we can ignore the change in effective pressure; therefore, the reconstructed seismic images reflect only those changes that arose from CO2 –water displacement. 3. Results 3.1. P-velocity and attenuation during the injection of CO2 in different states

(12)

The difference tomography can be applied to any kind of body wave. In the present study, we applied it to the first-arrived compressional wave (P-wave) obtained in the laboratory during a CO2 -injection test. Triangle mesh was applied to represent the velocity (as well as attenuation) field and to calculate ray paths. Velocity can be assigned to either every triangle cell (cell-mode) or every node (node-mode) of the mesh. The cell mode is useful to introduce discontinuity since the velocity in a cell is constant and the boundary between two neighbour cells may be discontinuous. Under the node mode, within each cell the velocity varies linearly with position, i.e., the velocity gradient is constant within each cell. In this paper, we selected the node mode and used a numerical two-dimensional wavefront propagator based on Huygens’ Principle to calculation ray path of the minimum travel time. Fig. 1(d) shows mesh geometry and ray paths calculated based on the wavefront migration method. In order to assure the reliability of the velocity models derived from the seismic tomography, a checkerboard model having 4x2 cells and 5% velocity anomalies was tested. The checkerboard model and the recovered image are shown in Fig. 1(e)–(f). The velocity anomalies were well imaged through the entire section with the exception of two end areas.

Fig. 3 shows an example of the waveforms and related basic data (travel times and wave amplitudes) corresponding to ray path S8–R8 during the injection of gaseous CO2 into the water-saturated sample. The arrival times have been converted to average ray-path velocity. This type of original data reveals a significant reduction in velocity and amplitude with increasing elapsed time. Figs. 4 and 5 show examples of tomography images for relative velocity and attenuation at selected times obtained during the injection of gaseous and liquid CO2 , respectively. The top and

2.4. Data-processing procedure First, data obtained under water-saturated conditions were selected as standard references for successive measurements during every injection of CO2 . The velocity distribution for every reference was determined on the basis of manually selected travel times. The absolute attenuation coefficient is not required for difference attenuation tomography; however, in order to convert the relative attenuation coefficient to Q for modelling purposes, we determined a mean value of Q in relation to a reference sample made of aluminium (with Q > 105 ) and having the same geometry

Fig. 3. Example of waveforms and basic data, including amplitude (Ap) and ray velocity (Vp) converted from travel times, corresponding to the horizontal path S8–R8 during the injection of gaseous CO2 into a water-saturated sample.

X. Lei, Z. Xue / Physics of the Earth and Planetary Interiors 176 (2009) 224–234

229

Fig. 4. Example of difference tomography images of relative velocity (dV) and the attenuation coefficient (dA) obtained during the injection of gaseous CO2 . The velocity distribution in the water-saturated sample is also shown. The labels below each pair of images indicate the elapsed time in minutes; e.g., the label “12-9” indicates that the associated images are relative velocity and the attenuation coefficient measured at 12 min relative to values measured at 9 min. “0” indicates the initial water-saturated condition (time = 0). In each pair, the left and right images correspond to P-velocity and P-attenuation, respectively. The overlapped thick lines in the right column indicate the well resolved area estimated from the checkerboard test.

Fig. 5. Difference tomography images of relative velocity (dV) and the attenuation coefficient (dA) obtained during the injection of liquid CO2 . The labels below each pair of images indicate elapsed time in minutes. See the caption to Fig. 4 for details.

230

X. Lei, Z. Xue / Physics of the Earth and Planetary Interiors 176 (2009) 224–234

Fig. 6. Temporal variation of the relative velocity (dV) at grids as indicated in Fig. 1(c). “A” and “B” correspond to high porosity regions at the lower and upper parts of the sample. “C” corresponds to a low porosity region at the middle part of the sample.

middle rows represent images of relative velocity and attenuation between two consecutive elapsed times. The bottom row represents images of relative velocity and attenuation at the given times with respect to the initial water-saturated condition (elapsed time = 0), i.e., the cumulative changes of velocity and attenuation. In all cases, upward CO2 migration was clearly imaged. We found that the spatial patterns of velocity reduction due to the replacement of water by CO2 are all similar. This feature is also true for increasing attenuation. By comparing these tomography images with the porosity distribution shown in Fig. 1, it is clearly demonstrated that regions of higher porosity resulted in a greater reduction in velocity and larger increase in attenuation. Fig. 6 shows temporal variation of relative velocities at three typical positions in the sample as indicated in Fig. 1(c) as “A”, “B” and “C”. “A” and “B” are positions of relatively higher porosity, while “C” corresponds to relatively lower porosity. P-velocity at “A” and “B” dropped with significantly larger amplitude. Fig. 7 shows that the final reductions in velocity are positively correlated with the local porosity. The distribution of porosity within the sample was estimated from the distribution of velocity under water-saturated condition based on a numerical correlation between P-velocity and porosity. Fig. 8 shows the mean reductions in velocity and the mean increase in attenuation coefficient with elapsed time, while Fig. 9 shows the mean reductions in velocity versus normalized increases in Q−1 . These mean values were calculated from the path-lengthweighted average of the relative tomography results. As shown in Fig. 8, the P-velocity initially decreased significantly before largely stabilizing. Attenuation showed an initial significant increase before largely stabilizing. During the injection of gaseous, liquid, and supercritical CO2 , the final velocity reductions were 7.5, 12, and 14.5%, respectively, while the attenuation coefficients increased by 20, 18, and 23 m−1 (Table 1). In converting the attenuation coeffi-

Fig. 7. Plot shows velocity reduction as a function of local porosity. The local porosity was estimated from the distribution of velocity in the water-saturated sample based on a linear correlation between P-velocity and porosity.

cient to Q−1 , we found that Q−1 for gas, liquid, supercritical states increased by factors of 3.3, 2.7, and 3.7 relative to water-saturated conditions. In summary, our results demonstrate that both the velocity and attenuation of P-waves are sensitive to partial saturation introduced by the displacement between the injected CO2 gas and pre-existing pore water; however, the change in attenuation of compressional waves is an order of magnitude greater than the change in velocity. 3.2. White and Dutta–Odé model As stated above, our experimental data demonstrate a large increase in attenuation owing to the replacement of water by CO2 fluid. Given that viscous loss is considered to be a major cause of this phenomenon, it is valuable to examine the WDO theory, as outlined in the Introduction, in the context of our experimental data. According to the WDO theory, the viscous loss in a partial saturation regime is modelled by considering porous rocks saturated with a single fluid (water in our case) but containing spherical regions filled with a second fluid (CO2 in our case). The WDO theory assumes that (1) the frequency-dependent stress field applied to the outer boundary of a unit cell is the same everywhere; i.e., patches (spheres) of heterogeneous saturation are much larger than the grain scale but much smaller than the wavelength, and (2) the basic cells do not interact. For mathematical simplicity, the idealized geometry of a unit cell consists of a gas-filled sphere of radius a placed at the centre of a water-saturated spherical shell of outer radius b (b > a). Minor overlap between neighbouring cells is introduced to

Fig. 8. Normalized average reductions in velocity and increases in attenuation versus elapsed time during the injection of gaseous, liquid, and supercritical CO2 . Note that before every injection of CO2 , imbibition was performed for more than 2 days to ensure complete saturation with water. Also shown are the fluid pressures at the inlet (Pin ) and outlet (Pout ), and confining pressure (Pc ).

X. Lei, Z. Xue / Physics of the Earth and Planetary Interiors 176 (2009) 224–234

231

Fig. 9. Change in Q−1 versus change in velocity caused by the injection of gaseous, liquid, and supercritical CO2 into a water-saturated sample. See text for details.

ensure that the total volume of the cells is equal to the rock volume. Thus, the gas saturation Sg is equal to (a/b)3 . This model describes wave velocity and attenuation as a function of gas saturation, wave frequency, and the distance between gas pockets, which can be considered as patch size. A full set of the formulas used in calculating the complex bulk modulus for a partially saturated porous rock can be found in the handbook by Mavko et al. (2003). Having described the basic concept behind the WDO theory, let us now turn to a comparison of the theory with our experimental data. First of all, water saturation must be estimated. We currently have no means by which to make direct measurements of the local gas saturation within the rock sample; however, we can estimate the average saturation either from the volume of water that flowed out from the sample or from the volume of CO2 that injected into the sample. This method is valid only for the period before the breakthrough (CO2 front migrated to the top boundary) was achieved. After the breakthrough, it is difficulty to make a precise estimation of gas saturation due to the fact that the volume fraction of CO2 and water flow out through the outlet was unknown. Second, it is better to convert the relative attenuation coefficient to absolute Q−1 . As mentioned in subsection 2.4, the average value of Q is 10 for the water-saturated sample. The relative attenuation coefficient dA can be then converted to dQ−1 using Eq. (12). Finally, we have to determine the material properties of the rock, water, and CO2 . All of the parameters necessary for the simulations are summarized in Table 2. The bulk modulus of the solid grains of Tako sandstone was chosen as an average of the bulk modulus of quartz and feldspar;

Fig. 10. Average strains against hydrostatic confining pressure during dry compression.

Ks = 38 GPa (Wang, 2000). The bulk and shear moduli of the solid matrix are estimated from seismic data obtained during dry compression. Given that the sample is characterized as a transverse isotropic medium and our results are associated with the minimum velocity Vp(0), propagating along the normal direction of the bedding plane. We used the bulk and shear moduli associated with Vp(0); Km = Gm = 6 GPa. The resulted Poison’s ratio of 0.125 is a reasonable value along the normal direction of the bedding plane. As is known, porosity decreases with increasing effective pressure, thus it is important to consider the compression effect. The correlation between porosity an effective pressure was estimated from the average volumetric strain data obtained during dry compression (Fig. 10). Since the effective pressure during all injections is less than 5 MPa, the reduction of porosity due to confining pressure is normally less than 0.3%. We calculated the permeability of the sample from our data. The resulting value is 15 mD, higher than the measured value of 7 mD reported by Xue and Ohsumi (2004a). The wave frequency was given as the dominant frequency 260 kHz obtained from the spectrum of waveform at partial saturation as shown in Fig. 3. After all of the material properties have been specified, the only free parameter is the patch size. We carried out numerical simulations with different b values and found that models with b = 1.3–1.5 mm yield results consistent with the observed velocity and attenuation. Fig. 11 shows the average velocity and attenua-

Table 2 Physical properties of Tako sandstone and fluids used in this study. Rock Ks Km Gm s ϕ

Bulk modulus of solid grains Bulk modulus of the solid matrix Shear modulus of the solid matrix Density of solid grains Average porosity Permeability

38 GPa 6 GPa 6 GPa 2500 g/m3 0.25 1.5 × 10−14 m2 (15 mD)

Fluid

Bulk modulus (GPa)

Density (kg/m3 )

Viscosity (Pa s)

Air Water Gaseous CO2 (23 ◦ C, 5 MPa) Liquid CO2 (23 ◦ C, 12 MPa) Supercritical CO2 (34 ◦ C, 12 MPa)

0.05 2.25 0.025a 0.066a 0.046a

1.29 997.0 124a 700a 620a

1.5e−5 1e−3 (23 ◦ C), 8e−4 (34 ◦ C) 1.7e−5 8.5e−5 7.0e−5

a

Xue and Ohsumi (2004b).

232

X. Lei, Z. Xue / Physics of the Earth and Planetary Interiors 176 (2009) 224–234

tion versus saturation, along with the theoretical values obtained using the WDO model for every CO2 injection. The well-known upper and lower bounds obtained using Gassmann’s relations are also plotted. To obtain the upper and lower bound, we used Voigt’s average and Reuss’s average to calculate the effective bulk modulus of the pore fluid (CO2 and water) in the Gassmann’s relations, respectively. Data after the breakthrough were also plotted in open circles but the saturation associated with these data was not estimated. The numerical results demonstrated that following an increase in CO2 saturation, attenuation of the compressional wave is characterized by a rapid increase, a peak at 30–40% CO2 saturation, and a gradual decrease until full CO2 saturation is achieved (if it is possible). Fig. 8 shows that at the final point, velocity was not completely

stabilized but maintained a gentle decreasing trend; the attenuation maintained a gentle increasing trend. Bearing this in mind as we return to Fig. 11, we believe that our CO2 -injection tests almost achieved the level corresponding to the peak attenuation as predicted by the WDO model. Thus, we conclude that the final obtained CO2 saturation is about 30–40% for all cases. There is no evidence of contrasting replace rates among the different states of CO2 . It is difficult to achieve a higher degree of CO2 saturation under the present experimental conditions. The estimated final saturation is in good agreement with the results calculated from X-ray CT data obtained during CO2 -injection tests on the same sandstone under similar injection conditions to those used in the present experiment.

Fig. 11. Comparison of observed velocity reductions and Q–1 with numerical results obtained using the White and Dutta–Odé model. The properties of the rock and fluid are described in Table 2. Open circles shows data obtained after the breakthrough (CO2 front migrated to the top boundary), the associated CO2 saturation was not estimated (see text for details).

X. Lei, Z. Xue / Physics of the Earth and Planetary Interiors 176 (2009) 224–234

233

4. Discussion It is verified that the difference tomography method, using differential arrival times and relative amplitudes of P-wave determined by cross-correlation, is efficient for imaging the replacing behaviour between injected CO2 and water in the pore space in rock sample in laboratory. It is possible to make reliable estimation of gas saturation with the use of time-lapsed tomography results of relative velocity and relative attenuation. Attenuation is especially important at gas saturations larger than 30–50% since P-velocity is not sensitive to gas saturation in high gas saturation regime. Laboratory-scale research under well-controlled conditions is helpful in understanding the mechanisms of fluid displacement between CO2 and water and in developing techniques applicable to field-scale problems; however, at field scales we have to use much lower frequencies. Direct measurements of the attenuation at fields are difficult to obtain in practice, and most commonly the low-frequency values must be extracted by extrapolating from measurements obtained at higher frequencies. For this reason, theoretical models play an important role. A good model should be valid over a wide frequency band. Our experimental results are in good agreement with results obtained using the White and Dutta–Odé model for partial saturation, indicating that mesoscopic fluid flow plays an important role in velocity dispersion and energy loss at ultrasonic frequencies. It is appropriate to mention the cut-off frequencies beyond which the WDO theory may become questionable (Dutta & Odé, 1979b). The first cut-off frequency (fI ) comes from the assumption that the wavelength in the medium is long compared to the dimension of the unit cell. The estimated b values of 1.3–1.5 mm are significant smaller than the wavelength (∼10 mm). However, in Dutta and Odé (1979b) they arbitrarily suggested that fI should be such that wavelength is at least 10 times larger than the intercell distance, which is 2b. For a b of 1.5 mm, the cut-off wavelength is 10 × 2 × 1.5 mm = 30 mm, above the wavelength of 10 mm. The second cut-off frequency (fII ) is inherent in Biot’s theory since the assumption of Darcy flow can no longer be justified for frequencies larger than fII given by the following relation (Biot, 1962), fII ≈

1 fc , 8

fc =

 2 fl

(13)

where fc is the transition frequency arose from the Biot’s theory, = porosity, fl = fluid density,  = viscosity of the pore fluid, = absolute permeability of the rock. By putting the material properties into Eq. (13), we can estimate the transition frequencies for different conditions. For water-saturated condition, the resulted transition frequency fc is 2.6 MHz and the corresponding cut-off frequency fII is 325 kHz – a value right in the range of frequencies used in this paper. For partial saturation a transition frequency would be slightly smaller than 2.6 MHz. As a result, the range of frequencies used in this paper is close to the cut-off frequencies of the WDO model. First of all, we should note that both fI and fII are arbitrarily chosen in Dutta and Odé (1979b) and could be somewhat conservative. Secondly, it is also important to note that the regular patchy geometry in WDO theory is idealized but helpful in identifying the effects of mesoscopic fluid flow. More realistic modelling should consider the pore-scale distribution in patchy saturation regime. The b values obtained by fitting experimental data with the WDO model can be considered as a kind of “mean” distance of patchy size distribution. Fig. 12 shows theoretical velocity reduction as a function of gas saturation for different b values. It is clearly demonstrated that for b values smaller enough the prediction of the WDO model is coincident with the lower bound of Gasmann’s equation for homogeneous saturation. While, increasing b value results in curves closing to the upper bound for patchy saturation. It is suggested that ultrasonic

Fig. 12. Comparison of velocity reductions obtained using the White and Dutta–Odé model with different patch distance (b). It is clearly demonstrated that for small patch distances the WDO model coincides with the lower bound from Gasmann’s equation for homogeneous saturation, while, greater b value gives a result closing to the upper bound from Gasmann’s estimation for patchy saturation.

laboratory conditions will generally not be described well with Gasmann’s equation (Mavko et al., 2003). Up to now, little experimental work has been performed for investigating the seismic properties of reservoir rocks saturated with CO2 and water. Our experimental data may shed light on the examination of applicability of existing models in ultrasonic frequencies. It is now appropriate to discuss the scaling problem associated with the WDO theory. As described above, for a given set of material properties of the rock and pore fluids, the velocity and attenuation of P-waves in the WDO model are described as a function of frequency and b at a given gas saturation (Sg). It has been reported that the frequency can apparently be treated as a geometry parameter in the form fb2 (Dutta & Odé, 1979a, 1979b). This scaling is exact when the model is treated in the manner of the original White model; however, it is only approximately true in the exact treatment of the corrected WDO model. In Fig. 13, we plotted the computed phase velocity and Q–1 as a function of fb2 for several gas saturations. The results calculated with fixed frequency (varying b) and fixed b (varying frequency) are plotted, but they overlap with each other. Most of the velocity dispersion, or the peak attenuation, is seen to occur for fb2 between 200 and 4000 for Sg = 0.9. For our experimental case, we have f = 260,000 Hz and b = 0.15 cm; thus, fb2 = 4400 Hz cm2 . Any combination of f and b that satisfies fb2 = 4400 Hz cm2 will result in the same velocity and attenuation as that obtained in the experiment. For a seismic wave with f = 100 Hz, b = 6.6 cm is required to observe the same results as those obtained in our experiment. If b is scale-independent, for f = 100 Hz we have fb2 is 22.5; in this case, reduction in velocity and increase in attenuation would be very small compared with our ultrasonic data. However, as is known the patch saturation or fingering is considered as a natural result due to capillary force-governed permeability and viscous diffusing (Homsy, 1987). Therefore, it is a very interesting issue to clearly describe the scale-dependence of b, i.e., patchy size. We believe that patchy saturation might show some kind of patchy size distribution, as many geological systems, such as cracks and faults, show a fractal size distribution over a wide range of sizes; the patchy size distribution may also be fractal. The injection of CO2 into a geological formation might result in the formation of CO2 pockets of widely varying size. If this assumption is true for real formations, the theoretical model derived from experimental data could be expanded easily to seismic frequency band. Attenuation in Massilon sandstone at the sonic frequency band (Murphy, 1982) can

234

X. Lei, Z. Xue / Physics of the Earth and Planetary Interiors 176 (2009) 224–234

however, the final obtained CO2 saturation is only about 30–40%. It is difficulty to achieve a higher degree of CO2 saturation. Acknowledgment This work was supported in part by Ministry of Economy, Trade and Industrial of Japan under the contract for “Research and Development of Underground Storage for Carbon Dioxide”. References

Fig. 13. Attenuation and velocity dispersion of compressional waves calculated using the White and Dutta–Odé model. The properties of the rock and fluid are described in Table 1; gaseous CO2 is considered.

be reproduced using the WDO model with a patch size of 5 cm. Additional studies, including laboratory experiments and field-based injection tests, are expected to further clarify the scale-dependence of patch saturation. It remains an interesting and important challenge to incorporate varying path sizes into the WDO theory. Such issues are key points of ongoing works. 5. Conclusions We carried out seismic measurements in the ultrasonic band during CO2 -injection experiments in porous sandstone samples under laboratory conditions. The obtained high-quality seismic data enabled detailed determination of the relative velocity and attenuation coefficient of compressional waves using the difference seismic tomography method. We confirmed that the difference tomography method, when employed using differential arrival times and relative amplitudes determined from the crosscorrelation of waveforms, is efficient in imaging the changes in velocity and attenuation associated with the partial replacement of water by CO2 . On average, as a result of gaseous, liquid, and supercritical CO2 partially replacing pore water, P-velocity fell by 7.5, 12, and 14.5%, respectively, while Q−1 increased by factors of 3.3, 2.7, and 3.7. Our results are in good agreement with the White and Dutta–Odé theory for partial saturation, indicating that fluid diffusion plays a major role in velocity dispersion and energy loss at ultrasonic frequencies. Numerical results obtained for patch sizes of b = 1.3–1.5 mm are consistent with the observed velocity and attenuation data. These results demonstrate that following an increase in CO2 saturation, compressional wave attenuation is characterized by a rapid increase, a peak at 30–40% CO2 saturation, and a gradual decrease until full CO2 saturation is achieved (if it is possible);

Aki, K., Richards, P.G., 2002. Quantitative Seismology, second edition. University Science Books, Sausalito, California, 700 pp. Arts, R., Eiken, O., Chadwick, A., Zweigel, P., van der Meer, L., Kirby, G., 2004. Seismic monitoring at Sleipner underground CO2 storage site (North Sea). In: Baines, S.J., Worden, R.H. (Eds.), Geological Storage of Carbon Dioxide, 233. Geological Society, London, pp. 181–191 (special publication). Benson, S.M., 2005. Overview of geological storage of CO2 . In: Thomas, C., Benson, S.M. (Eds.), Carbon Dioxide Capture for Storage in Deep Geologic Formations, vol. 2, pp. 665–672. Biot, M.A., 1962. Mechanics of deformation and acoustic propagation in porous media. J. Appl. Phys. 33, 1482–1498. Carcione, J.M., Picotti, S., Gei, D., Rossi, G., 2006. Physics and seismic modeling for monitoring CO2 storage. Pure Appl. Geophys. 163, 175–207. Chadwick, R.A., Eiken, O., Lindeberg, E., 2004. 4D geophysical monitoring of the CO2 plume at Sleipner, North Sea: Aspects of uncertainty. In: Proceedings of the 7th SEGJ International Symposium, Sendai, pp. 24–26. Dutta, N.C., Odé, H., 1979a. Attenuation and dispersion of compressional waves in fluid-filled porous rocks with partial gas saturation (White model) – Part I: Biot theory. Geophysics 44, 1777–1788. Dutta, N.C., Odé, H., 1979b. Attenuation and dispersion of compressional waves in fluid-filled porous rocks with partial gas saturation (White model) – Part II: Results. Geophysics 44, 1789–1805. Dutta, N.C., Odé, H., Seriff, A.J., 1979. On White’s model of attenuation in rocks with partial gas saturation. Geophysics 44, 1806–1812. Harris, J.M., Nolen-Hoeksema, R.C., Van Schaack, M., Lazaratos, S.K., Rector, J.W., 1995. High-resolution crosswell imaging of a west Texas carbonate reservoir: Part I – Project summary and interpretation. Geophysics 60, 667–681. Homsy, G.M., 1987. Viscous fingering in porous media. Annual review. Fluid Mechanics 19, 271. International Panel on Climate Change (IPCC), Climate change 1995. The Science of Climate Change, Cambridge University Press, Cambridge, UK, 584 pp. International Panel on Climate Change (IPCC), Climate Change 2001. The Scientific Basis: Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change (Climate Change 2001), Cambridge University Press, Cambridge, UK, 892 pp. Li, G., 2003. 4D seismic monitoring of CO2 flood in thin fractured carbonate reservoir. The Leading Edge, 690–695. Lin, W., Nakamura, T., Fujita, M., 1997. Variation of ultrasonic wave amplitude in Inada granite and Tako sandstone under triaxial compression process:. In: 32nd Japanese National Conference on Geotechnical Engineering, pp. 1249–1250. Lo, W.-C., Sposito, G., 2005. Wave propagation through elastic porous media containing two immiscible fluids. Water Resources Res. 41, W02025, doi:10.1029/2004WR003162. Lu, J.F., Hanyga, A., 2005. Linear dynamic model for porous media saturated by two immiscible fluids. Int. J. Solids Struct. 42, 2689–2709. Murphy, W.F., 1982. Effect of partial water saturation on attenuation in Massilon sandstone and Vycor porous glass. J. Acoust. Soc. Am. 71, 1458–1468. Mavko, G., Mukerji, T., Dvorkin, J., 2003. The Rock Physics Handbook. Cambridge University Press, Cambridge, UK, 329 pp. Saito, H., Azuma, H., Nobuoka, D., Tanase, D., Xue, Z., 2006. Time-lapse crosswell seismic tomography monitoring of CO2 at an onshore aquifer, Nagaoka, Japan. Butsuri-Tansa 59, 30–36. Thomsen, L., 1986. Weak elastic anisotropy. Geophysics 51, 1954–1966. Tuncay, K., Corapcioglu, M.Y., 1996. Body waves in poroelastic media saturated by two immiscible fluids. J. Geophys. Res. 101 (25), 149–159. Wang, Z., 2000. The Gassmann equation revisited: comparing laboratory data with Gassmann’s predictions. In: Wang, Nur (Eds.), Seismic and Acoustic Velocities in Reservoir Rocks, vol. 3, Recent Developments, pp. 8–23. Walsh, J.B., 1995. Seismic attenuation in partially saturated rock. J. Geophys. Res. 100, 15,407–15,424. White, J.E., 1975. Computed seismic speeds and attenuation in rocks with partial gas saturation. Geophysics 40, 224–232. Xue, Z., Lei, X.-L., 2006. Laboratory study of CO2 migration in water-saturated anisotropic sandstone, based on P-wave velocity imaging. Exploration Geophys./Butsri-Tansa/Mulli-Tamsa 59, 10–18. Xue, Z., Ohsumi, T., 2004a. Laboratory measurements on gas permeability and Pwave velocity in two porous sandstones during CO2 flooding. Shigen-to-Sozai 120, 91–98 (in Japanese). Xue, Z., Ohsumi, T., 2004b. Seismic wave monitoring of CO2 migration in watersaturated porous sandstone. Exploration Geophys. 35, 25–32.