© L. Li & J. Ma, published by IFP Energies nouvelles, 2019

1 Introduction

The 4D seismic method is a remarkable approach with which to monitor the distribution of CO₂ plume in the reservoir during CO₂-EOR or aquifer storage (Ivanova et al., 2012; White, 2009, 2013). The map of CO₂ plume from 4D seismic monitoring data is used to evaluate the dead oil zone and the efficiency of CO₂ flooding. The estimation of CO₂ volume in the reservoir is our ultimate goal and one of the most powerful evidence to prove the safety of CO₂ geological storage.

Conventional 4D seismic interpretation methods (Avseth et al., 2006) assume that porosity never be changed before and after fluids (CO₂ and water) injection. The bulk modulus of the mineral matrix and bulk modulus of the porous rock frame will not be changed as well. Then fluids substitution methods based on Gassmann’s equation is used to make rock physics model (Mavko et al., 1998), and further to make 4D seismic forward models for interpretation.

However, more scientific evidence proved that reservoir pore system or pore structure will be changed after CO₂ injection. Berrezueta et al. (2013) observed that intergranular clay matrix detachment and partial removal in CO₂ flooded sandstone samples that lead to porosity increased 3.75%. Geochemical data observed in CO₂-injected fields (Hovorka, 2009; Karamalidis et al., 2010; Kharaka et al., 2006; Thordsen et al., 2012) indicate that injection lowers the pH of the brine by dissolving a few tens of milligrams per litre of calcium (Ca), iron (Fe), and manganese (Mn) from the matrix as the CO₂ front moves through the rock-water system; this phenomenon also indicates that dissolution occurs not only in the injection well but also in the reservoir. Tambach et al. (2015) simulated CO₂ dissolution in the K12-B area in Holland and found that CO₂ dissolution can increase the porosity and decrease the formation pressure. Vanorio (2015) performed time-lapse petrophysical experiments to observe the effects of fluid–rock interactions on petrophysical parameters and observed that the P- and S-wave velocities of rocks decreased before and after CO₂ injection under the same pressure regardless of whether the injection was performed in a carbonate or a sandstone; this phenomenon is a result of the coupling between compaction and variations in the porosity. Fӧster et al. (2010) observed that approximately 1–5% of the porosity in the Ketzin CO₂ sequestration area in Germany resulted from CO₂ dissolution of clastic rock. Lu et al. (2016) experimentally examined the dissolution effect of CO₂ on deposits and found that CO₂ dissolution could increase the porosity and that the magnitude of the increase was associated with the pore throat radius. Several authors also found the porosity decreased after CO₂ flooding (Skorpa et al., 2017).

Establishing a reasonable 4D seismic forward model is the key to correctly interpreting 4D seismic data, and 4D seismic responses result from variations in the fluid saturation or pressure in a deposit. Accordingly, the goals of using the 4D seismic method in this manner are to identify the effects of the CO₂ saturation and pressure, to monitor the state of underground CO₂, and to provide a theoretical understanding of the safety of geological CO₂ sequestration.

The formation pressure and fluid saturation can vary after CO₂ injection; in addition, the porosity also changes due to the fracturing either induced by high-pressure CO₂ injection or due to CO₂ dissolution of the rock. Thus, it is necessary to establish a 4D seismic forward model that considers variations in the fluid saturation, porosity, and formation pressure simultaneously, and accurate estimations of P- and S-wave velocities constitute the basis of establishing a 4D seismic forward model.

For present-day 4D seismic exploration, many P- and S-wave velocities estimation models are based on the Hertz-Mindlin model (Duffaut and Landrø, 2007; Mindlin, 1949; Saul and Lumley, 2015). This method considers the effects of both pressure and fluid. Duffaut and Landrø (2007) proposed the use of the Hertz-Mindlin expression to calculate P- and S-wave velocities and establish 4D seismic forward and inverse models in poorly consolidated sandstone. Saul and Lumley (2015) developed an inelastic model based on diagenesis and grain cement using the Hertz-Mindlin expression to describe the pressure sensitivity of rock.

However, none of these methods fully considers the effects of fluids. Instead, we propose that the effects of fluids should be considered in two aspects: (1) the effects of variations in the fluid composition, which can be achieved by the substitution of different fluids in the Gassmann equation, and (2) the effects of variations in the porosity on the seismic response after CO₂ injection, which causes changes in the fluid. Ignoring variations in the porosity could result in misinterpretation; for example, the main factor influencing the seismic response could be the fluid composition, but neglecting porosity variations could result in the conclusion that pressure is the main factor.

The Digby equation assumes that cemented rock grains are present, whereas the Hertz-Mindlin expression is applicable to uncompacted sandstone reservoirs. The strata of interest in our study area, the Weyburn oil reservoir, can be divided into a marly layer and a vuggy layer; the former is composed of granular limestone, while the latter is finely crystalline to granular dolomite. Based on a petrological analysis and Scanning Electron Microscopy (SEM) inspection (Jensen, 2016; Njiekak et al., 2013), the Digby equation is suitable in the study area. Here, we improve upon the Digby equation to calculate the bulk modulus and shear modulus of dry rock and to estimate the P- and S-wave velocities of saturated rock under different pressures in combination with the Gassmann equation. The P- and S-wave velocities both pre- and post-CO₂ injection are the basis of a 4D seismic forward model. The S-wave velocity is sensitive to the pressure rather than the fluid composition, whereas the P-wave velocity is affected by both the pressure and the fluid composition. Therefore, we focus mainly on variations in the P-wave velocity in this study. A 4D seismic forward model is established based on this relationship to investigate the effects of changes in the fluid composition and pressure on the Amplitude Versus Offset (AVO) gradient and intercept. In addition, the CO₂ conditions studied here can be applied not only in the Weyburn deposit but also in other areas as well, and our motivation is to provide new 4D seismic interpretations.

2 Geological and geophysical background

The Weyburn deposit was discovered in 1954, and water was injected in 1965; more recently, CO₂-Enhanced Oil Recovery (EOR) was initiated in 2000. The deposit is located in the Midale layer at depths of 1300–1500 m. The Midale layer can be divided into two layers: an upper marly dolomite (2–12 m) and a lower vuggy limestone (8–22 m). The marly layer has a high porosity (29%) and a low permeability (mean (M) = ~10 mD), whereas the vuggy layer has a lower porosity (10%) but a higher permeability (M = ~50 mD) (Brown, 2002). When injecting CO₂, the injection pressure at the bottom of the well is 22 MPa, the original pore pressure is 15 MPa, and the pressure at the bottom of the production well is controlled to be 7–8 MPa (Ma and Morozov, 2010; White, 2013). According to research on the Weyburn deposit conducted by the British Geological Survey (Riding and Rochelle, 2005), the high-pressure injection of CO₂ can result in fracturing and increase the total porosity by 0.5–1.9%.

Well 4-23 was drilled at the edge of Phase 1A in August 2000, approximately one month before CO₂ injection (Fig. 1), to obtain well log curves (Fig. 2). The S-wave velocity data originate only from this well, which is therefore treated as the basis of the AVO model and deposit estimates. However, cores were collected from Wells 14-11 and 14-23 rather than from Well 4-23. Ultrasonic measurements were performed on four core samples from the marly and vuggy units (Brown, 2002). Their density, porosity, and permeability values are shown in Table 1. Our study is also based on these results.

3 Theory

Digby (1981) proposed a method to calculate the bulk modulus and shear modulus of closely compacted dry rock based on changes in the pressure:

(1)

(2)where b can be expressed as:

(3)while d satisfies equation (4):

(4)where K _dry and μ _dry are the bulk and shear moduli, respectively, of dry rock; υ and μ _ma are Poisson’s ratio and the shear modulus, respectively, of the solid grains; Φ is the porosity; C _p is the coordination number; p is the differential pressure; a is the radius of the small, flat, circular regions of neighbouring bonded particles before deformation; b is the radius after deformation; and R is the grain radius. K _dry, μ _dry and μ _ma have the same unit of measurement (GPa), while the unit of measurement for p is MPa, and υ, φ, C _p , a/R and b/R are all dimensionless. The bulk and shear moduli of the rock matrix are calculated by using the Voigt-Reuss-Hill (V-R-H) method.

Murphy (1982) assumed C _p to be a constant, summarized C _p for different porosities, and concluded that C _p is proportional to e ^1−ϕ . The constant C _p values for different porosities summarized by Murphy are shown in Table 2 (where the porosity and C _p are dimensionless).

Table 2 shows the C _p values corresponding to the porosities summarized by Murphy, where the porosity is listed discretely as 0.25, 0.3, 0.35,…, 0.7. In Figure 3, which displays a fitting diagram based on Murphy’s statement that C _p is proportional to e ^1−ϕ , we compare C _p with e ^1−ϕ and use the values obtained from linear fitting for subsequent calculations using the equation C _p  = 11.759e ^1−ϕ −12.748.

After fitting data with this formula, the corresponding C _p values for a continuous change in the porosity can be calculated. However, C _p varies with pressure (Wang et al., 2018). Based on this linear relationship, we propose a revised equation (5) for C _p :

(5)where W is a weighting factor that represents several effects, such as the effects of the pressure and lithology. The P-wave velocity is usually measured in actual oil fields, and the best W value is obtained by subtracting the absolute value of the P-wave velocity in the actual logging data from the predicted P-wave velocity and establishing a minimum or zero difference (Eq. (6)):

(6)where is a measured P-wave velocity and (W) is a predicted value; the units of both and (W) are km/s. Substituting equation (5) into equations (1) and (2), we can obtain K _dry(W) and μ _dry(W), which can then be substituted into Gassmann’s equation, following which we can obtain the bulk and shear moduli of a saturated rock that contains W. Then, (W) can be represented as equation (7), where K _sat(W) and μ _sat(W) are the bulk and shear moduli, respectively, of a saturated rock that contains W; the units of both K _sat(W) and μ _sat(W) are GPa, and ρ is the density expressed in units of g/cm³. Equation (8) is an expansion of equations (6) and (7). Finally, we can calculate the shear modulus μ _dry by substituting the predicted W and C _p values into equation (2) to predict the S-wave velocity:

(7)

(8)

4 Calculation of elastic wave velocities

4.1 Method validation

As early as 2002, Brown performed the petrophysical tests on the rock of Weyburn field’s core samples in Colorado School of Mines. Four core samples were acquired from Marly and Vuggy units in two wells. For these four core samples, Brown measured the wave velocity once when pressure increased from low to high until the maximum pressure was reached. Then the velocity was measured again when pressure decreased from high to low. Two groups of wave velocity under two different pressures were obtained (Brown, 2002). For the compressional wave velocity, the average of the two measurements was selected. For the shear-wave velocity, because two values of V _s1 and V _s2 were measured for each core sample and each value of V _s1 or V _s2 was measured two times, there were four shear-wave velocities for each core sample. We took their average as the shear wave velocity in an isotropic medium and compared it with the predicted values (Brown, 2002).

Variations in the P-wave velocity during geological CO₂ sequestration result from changes in both the fluid composition and the fluid pressure. To identify the influences of these two effects, the marly layer, which has the higher porosity between the two units, is selected as the subject of this study. Two marly rock cores, PC12-13 1410 and PC14-11 1420.2, are used to verify our method. Then, calculation models are developed to estimate the elastic wave velocities based on the rock testing results of rock core PC12-13 1410.

First, the S-wave velocity is estimated using the Digby equation in combination with the Gassmann expression based on the measured P-wave velocity, as shown in Tables 3 and 4 and Figures 4 and 5. The values of the weight coefficient W in Tables 3 and 4 are calculated using equation (8). The absolute errors of the estimated S-wave velocities are less than 5%, resulting in a prediction curve with a high accuracy, thereby indicating the efficiency of this method.

4.2 P- and S-wave velocity calculations

After verifying the accuracy of the proposed method, the physical parameters of rock core PC12-13 1410 are used for the following analysis because its average porosity (29%) is similar to that of the marly layer.

The P-wave velocity is calculated using the Digby equation with variable porosity and pressure. The Digby equation includes the pressure effect, and both the Digby equation and the coordination method proposed in our study include the porosity effect, which constitutes the theoretical basis of our study. Notably, based on the values of the weight coefficient shown in Tables 3 and 4, the coefficient remains constant under similar pressures with varied porosities (the porosities of these two rock cores are 0.29 and 0.33), implying that the weight coefficient is associated only with the pressure when the lithology is the same. On the other hand, the coordination number is associated with both the pressure and the porosity, enabling the S- and P-wave velocities of rock core PC12-13 1410 to be calculated under different porosities.

The elastic parameters employed herein are from Ma and Morozov (2010). To infer the effects of variations in the fluid on the P-wave velocity, P-wave velocities are calculated under different pressures, different porosities, and different CO₂ saturations, as shown in Figure 6, the results of which reflect the variation in the P-wave velocity with a changing porosity. Different shapes represent different pore pressures, and different colours represent different CO₂ saturations. The results show that P-wave velocity at various pressures can be distinguished under the different porosity and CO₂ saturation. However, for a given pressure with variable porosity, the P-wave velocities under different CO₂ saturations can overlap. For example, when the pore pressure is 8 MPa and the porosity is less than 0.24, the P-wave velocity under 5% CO₂ saturation is greater than that under 20% CO₂ saturation; however, when the porosity is greater than 0.24, the velocity pattern is reversed. This relationship implies that the effect of the porosity must be considered when examining variations in the fluid saturation.

We focus on the effects of porosity variations on the S- and P-wave velocities and on the forward models before and after CO₂ injection, following which P-wave velocities under different porosities and saturations are calculated for two fixed pressures. As shown in Figure 7, as the porosity varies, the P-wave velocities under different pressures cannot be identified because they overlap one another, that is, as the porosity increases by 1%, the P-wave velocity for a pore pressure of 8 MPa and a CO₂ saturation of 10% (red triangle) overlaps the P-wave velocity for a pore pressure of 15 MPa and a CO₂ saturation of 5% with the original porosity (black star). For a pore pressure of 8 MPa with increasing porosities and CO₂ saturations, the P-wave velocities cannot be distinguished from those at a pore pressure of 15 MPa, highlighting the importance of considering porosity variations in the forward model. In Figure 6, the main factor influencing the P-wave velocity is still the fluid pressure; however, with varying porosity (Fig. 7), the P-wave velocities of two pressures cannot be distinguished from one another.

Figure 8 shows the P-wave velocities at pore pressures of 8 MPa and 15 MPa. Similar to Figure 7, when the porosity is increased by 1% (red triangles), the P-wave velocity for a pore pressure of 8 MPa overlaps that for a pore pressure of 15 MPa.

The P-wave velocity varies with the porosity. However, the effects of porosity variations on the S-wave velocity remain to be studied. Here, we study these effects with different pressures and CO₂ saturations. In Figure 9, the x-axis represents the porosity, and the y-axis is the S-wave velocity; the figure shows data for different pressures and saturations and porosities. Clearly, the variation in the S-wave velocity is determined by the fluid pressure.

Here, we study the changes in the S-wave velocity with changes in the porosity. In Figure 10, the variation in the S-wave velocity is determined mainly by the pressure, although S-wave velocities at different pressures can still overlap in response to large variations in the porosity at high CO₂ saturations (purple triangles). These results also demonstrate the significance of considering both the CO₂ saturation and the porosity in the forward model.

5 Establishing a 4D seismic forward model

Based on the S- and P-wave velocities obtained under different pressures, CO₂ saturations and porosities, the reflection coefficients are calculated using the Zoeppritz equation to compute the synthetic seismogram through the angle domain from convolution of the Ricker wavelet.

5.1 Two-layer model

A two-layer model is first developed for synthetic seismograms at 8 MPa with different saturations and porosities, and difference profiles are obtained by comparing these synthetic seismograms (Fig. 11) with the parameters shown in Table 5. As shown in Figure 11, the differences in the synthetic seismograms increase under similar pressures and saturations but with increasing porosity.

The AVO gradient intercept figure (Fig. 12) corresponds to the two-layer model obtained above. As shown in Figure 12, variation in the porosity leads to variation in the AVO intercept, and variation in the CO₂ saturation results in a change in the AVO gradient, which is consistent with the responses in the two-layer seismic model.

To investigate the effects of porosity variations on the seismic response and AVO, synthesize seismograms for two-layer models with different porosities and pressures are established (Fig. 13) with the parameters shown in Table 6. In Figure 13a, the AVO trend of the difference between the synthetic seismograms increases in amplitude with an increase in the angle, whereas in Figure 13d, such a change is not obvious, implying the importance of considering variations in the porosity.

To directly observe the effects of porosity and pressure variations on the AVO, the AVO gradient intercept is calculated based on the two-layer model data in Table 6. Figure 14 shows that variation in the fluid saturation leads to variation in the AVO gradient, where the variation in the AVO intercept is due to the variations in the pressure and porosity.

5.2 Well model

Based on the abovementioned two-layer models, a well model can be developed using data from Well 4-23 to compute synthetic seismograms and the differences between those seismograms. The well log curve from Well 4-23 is shown in Figure 2. First, a well model is derived using the original well log data to create a synthetic seismogram; then, the fluid is substituted based on the proposed method to predict the P- and S-wave velocities after CO₂ injection to fully develop the well model. Before fluid substitution, the pore pressure is 15 MPa, and the CO₂ saturation is 0; after fluid substitution, the pore pressure is 8 MPa, and the CO₂ saturation is 20%. To consider porosity variations, two synthetic seismograms are generated after fluid substitution: one without changes in the porosity and one in which the porosity is increased by 1.5%. Figures 15a–15c show the synthetic seismograms of the well model before and after fluid substitution. The differences between the synthetic seismograms with varying porosity are greater than those between the synthetic seismograms with a consistent porosity, further demonstrating the effect of a changing porosity on the seismic response during the establishment of the forward model. However, further discussion of the unclear difference in Figure 15c is required.

6 Discussion

6.1 The effects of fluid variations on the P-wave velocity

Figure 7 shows that, at a constant pressure, the P-wave velocity decreases and then increases with increasing CO₂ saturation and porosity. We propose two explanations for this observed pattern. First, the porosity increases to a critical value, at which point the carrier of the seismic wave converts from minerals to fluid (Mavko et al., 1998). When the porosity is lower than the critical value, the minerals impose a dominant influence on the P-wave velocity; therefore, with increasing CO₂ saturation, the P-wave velocity decreases. After reaching the critical porosity, the main influence is derived from the fluid, and the variation in the P-wave velocity is associated with the properties of this fluid. Second, the properties of the fluid mixture vary with pressure. The bulk moduli and densities of the fluid mixture at different pressures are shown in Figures 16a and 16b, respectively. After fluid substitution, the pore pressure of the deposit is 8 MPa, which represents the bubble point. In Figure 16, the bulk moduli and densities of the fluid mixture exhibit their greatest changes below 8 MPa in each state. Due to great variations in the fluid properties, the P-wave velocity can vary substantially at high porosity values.

6.2 Analysis of the synthetic seismogram of the well model

According to the synthetic seismogram of the well model presented in Section 5.2, a distinct difference is observed between the original seismic traces and the seismic traces after fluid substitution, implying that CO₂ injection causes variations in the fluid saturation and pore pressure, which can be identified based on the seismic response. However, no significant difference exists between the synthetic seismograms with and without considering porosity variations (Fig. 15c), and the porosity effect shown in the synthetic seismogram difference profile from the well model is not as clear as that from the two-layer model; this discrepancy is likely due to two reasons. First, according to the well log curve in Figure 2, the marly layer has a higher porosity than the vuggy layer. Since the seismic responses associated with porosity variations are due to the fluid effect, the porosity below 0.1 cannot significantly influence the P- and S-wave velocities or the seismic response, even if the porosity varies by 1.5%. Second, the deposit layer is thin, the marly layer is even thinner. Therefore, any variation in the marly layer is hardly observable due to the sampling restrictions of the preceding depth-time conversion procedure when synthesizing the seismogram.

7 Conclusion

Based on this study, we can conclude the following:

1. Through calculations of the P- and S-wave velocities of the marly layer in the Weyburn deposit under different conditions, variations in the porosity can result in variations in the P- and S-wave velocities and even cause the velocities to overlap at different pressures. Therefore, in addition to variations in the pressure and fluid saturation, changes in the porosity must be considered when estimating 4D seismic P- and S-wave velocities, especially when the porosity of the deposit is relatively high.

2. According to two layer synthetic seismogram, when comparing cases with and without porosity variations before and after fluid substitution, porosity variations may lead to greater seismic responses and can even change the AVO trend of the synthetic seismogram.

3. Variations in porosity can provide a new possible explanation for observed 4D seismic patterns and can reduce the diversity of interpretations.

A 4D seismic forward model has been established. Inversion of 4D seismic data and interpretation 4D seismic data will be done on the basis of this 4D seismic forward model, which provide basis for the safety of CO₂ geological storage.

References

All Tables

All Figures