Regular Article
The influence of pore system change during CO_{2} storage on 4D seismic interpretation
National & Local Joint Engineering Research Center of Carbon Capture and Storage Technology, Department of Geology, Northwest University, 710069 Xi’an, PR China
^{*} Corresponding author: lilin_hyg@163.com
Received:
23
April
2019
Accepted:
16
August
2019
A 4D seismic forward model constitutes the foundation of 4D seismic inversion. Here, in combination with the Gassmann equation, the Digby model is improved to calculate the Swave velocity, and the resulting equation is verified using rock testing results. Then, considering the influences of changes in the pore pressure, CO_{2} saturation and porosity on the P and S wave velocities, rock testing results from a CO_{2} injection area in the Weyburn field are used to calculate the P and Swave velocities of the reservoir. These P and Swave velocities are found to overlap under different pressure conditions with or without considering porosity variations. Therefore, twolayer models and well models are developed to simulate synthetic seismograms; the models considering porosity variations may provide greater seismic responses and different Amplitude Versus Offset (AVO) trends in the synthetic seismogram profiles than those without considering porosity variations. Thus, porosity variations must be considered when establishing 4D seismic forward models.
© L. Li & J. Ma, published by IFP Energies nouvelles, 2019
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The 4D seismic method is a remarkable approach with which to monitor the distribution of CO_{2} plume in the reservoir during CO_{2}EOR or aquifer storage (Ivanova et al., 2012; White, 2009, 2013). The map of CO_{2} plume from 4D seismic monitoring data is used to evaluate the dead oil zone and the efficiency of CO_{2} flooding. The estimation of CO_{2} volume in the reservoir is our ultimate goal and one of the most powerful evidence to prove the safety of CO_{2} geological storage.
Conventional 4D seismic interpretation methods (Avseth et al., 2006) assume that porosity never be changed before and after fluids (CO_{2} 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_{2} injection. Berrezueta et al. (2013) observed that intergranular clay matrix detachment and partial removal in CO_{2} flooded sandstone samples that lead to porosity increased 3.75%. Geochemical data observed in CO_{2}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_{2} front moves through the rockwater 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_{2} dissolution in the K12B area in Holland and found that CO_{2} dissolution can increase the porosity and decrease the formation pressure. Vanorio (2015) performed timelapse petrophysical experiments to observe the effects of fluid–rock interactions on petrophysical parameters and observed that the P and Swave velocities of rocks decreased before and after CO_{2} 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_{2} sequestration area in Germany resulted from CO_{2} dissolution of clastic rock. Lu et al. (2016) experimentally examined the dissolution effect of CO_{2} on deposits and found that CO_{2} 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_{2} 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_{2} saturation and pressure, to monitor the state of underground CO_{2}, and to provide a theoretical understanding of the safety of geological CO_{2} sequestration.
The formation pressure and fluid saturation can vary after CO_{2} injection; in addition, the porosity also changes due to the fracturing either induced by highpressure CO_{2} injection or due to CO_{2} 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 Swave velocities constitute the basis of establishing a 4D seismic forward model.
For presentday 4D seismic exploration, many P and Swave velocities estimation models are based on the HertzMindlin 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 HertzMindlin expression to calculate P and Swave 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 HertzMindlin 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_{2} 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 HertzMindlin 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 Swave velocities of saturated rock under different pressures in combination with the Gassmann equation. The P and Swave velocities both pre and postCO_{2} injection are the basis of a 4D seismic forward model. The Swave velocity is sensitive to the pressure rather than the fluid composition, whereas the Pwave velocity is affected by both the pressure and the fluid composition. Therefore, we focus mainly on variations in the Pwave 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_{2} 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_{2}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_{2}, 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 highpressure injection of CO_{2} can result in fracturing and increase the total porosity by 0.5–1.9%.
Well 423 was drilled at the edge of Phase 1A in August 2000, approximately one month before CO_{2} injection (Fig. 1), to obtain well log curves (Fig. 2). The Swave 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 1411 and 1423 rather than from Well 423. 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.
Fig. 1 The Phase 1A area of CO_{2} sequestration in the Weyburn field, Canada. Well 423 is the logging well, whereas Well 1213 and Well 1411 are the coring wells. 
Fig. 2 Well log data from Well 423. From left to right are the Pwave velocity (V _{ p }), Swave velocity (V _{ s }), resistivity (Rt) log, Gamma Ray (GR) log, shale content (V _{sh}), porosity (Por; the dashed line indicates the effective porosity, while the solid line indicates the total porosity), and density. 
Properties of the rock cores.
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:
(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 VoigtReussHill (VRH) 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).
Summary of information on various porosities and their corresponding C _{ p } values.
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.
Fig. 3 Linear fitting between C _{ p } and e ^{1−φ }. Black dots are Murphy’s summarization. Black line is result of linear regression. 
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 Pwave velocity is usually measured in actual oil fields, and the best W value is obtained by subtracting the absolute value of the Pwave velocity in the actual logging data from the predicted Pwave velocity and establishing a minimum or zero difference (Eq. (6)):
(6)where is a measured Pwave 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^{3}. 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 Swave velocity:
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 shearwave 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 shearwave 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 Pwave velocity during geological CO_{2} 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, PC1213 1410 and PC1411 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 PC1213 1410.
First, the Swave velocity is estimated using the Digby equation in combination with the Gassmann expression based on the measured Pwave 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 Swave velocities are less than 5%, resulting in a prediction curve with a high accuracy, thereby indicating the efficiency of this method.
Fig. 4 Comparison between the estimated and measured Swave velocities in PC1213 1410 from the marly layer. The error bar shows a 95% confidence interval for the mean of measured Swave velocity. 
Fig. 5 Comparison between the estimated and measured Swave velocities in PC1411 1420.2 from the marly layer. The error bar shows a 95% confidence interval for the mean of measured Swave velocity. 
Results for the marly layer in PC1213 1410.
Results for the marly layer in PC1411 1420.2.
4.2 P and Swave velocity calculations
After verifying the accuracy of the proposed method, the physical parameters of rock core PC1213 1410 are used for the following analysis because its average porosity (29%) is similar to that of the marly layer.
The Pwave 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 Pwave velocities of rock core PC1213 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 Pwave velocity, Pwave velocities are calculated under different pressures, different porosities, and different CO_{2} saturations, as shown in Figure 6, the results of which reflect the variation in the Pwave velocity with a changing porosity. Different shapes represent different pore pressures, and different colours represent different CO_{2} saturations. The results show that Pwave velocity at various pressures can be distinguished under the different porosity and CO_{2} saturation. However, for a given pressure with variable porosity, the Pwave velocities under different CO_{2} saturations can overlap. For example, when the pore pressure is 8 MPa and the porosity is less than 0.24, the Pwave velocity under 5% CO_{2} saturation is greater than that under 20% CO_{2} 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.
Fig. 6 Variations in the Pwave velocity with different porosities under different pressure and saturation conditions. The xaxis is the porosity, and the yaxis is the Pwave velocity. The triangle symbols represent a pore pressure of 8 MPa, and the star symbols represent a pore pressure of 15 MPa. The black, green, blue, and red colours represent the Pwave velocities associated with CO_{2} saturations of 5%, 10%, 15%, and 20%, respectively. 
We focus on the effects of porosity variations on the S and Pwave velocities and on the forward models before and after CO_{2} injection, following which Pwave velocities under different porosities and saturations are calculated for two fixed pressures. As shown in Figure 7, as the porosity varies, the Pwave velocities under different pressures cannot be identified because they overlap one another, that is, as the porosity increases by 1%, the Pwave velocity for a pore pressure of 8 MPa and a CO_{2} saturation of 10% (red triangle) overlaps the Pwave velocity for a pore pressure of 15 MPa and a CO_{2} saturation of 5% with the original porosity (black star). For a pore pressure of 8 MPa with increasing porosities and CO_{2} saturations, the Pwave 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 Pwave velocity is still the fluid pressure; however, with varying porosity (Fig. 7), the Pwave velocities of two pressures cannot be distinguished from one another.
Fig. 7 Variations in the Pwave velocity under different porosity, pressure and saturation conditions. The xaxis is the variation of porosity, and the yaxis is the Pwave velocity. The triangle symbols represent a pore pressure of 8 MPa, and the star symbols represent a pore pressure of 15 MPa. The black colour represents the Pwave velocity trend at a CO_{2} saturation of 5% and a porosity of 0.18–0.29, while the red, blue, green, grey and purple colours represent the Pwave velocity trends at CO_{2} saturations of 10%, 15%, 20%, 25%, and 30% and porosities that are 1–5%, respectively, higher than the data represented by the black colour. 
Figure 8 shows the Pwave velocities at pore pressures of 8 MPa and 15 MPa. Similar to Figure 7, when the porosity is increased by 1% (red triangles), the Pwave velocity for a pore pressure of 8 MPa overlaps that for a pore pressure of 15 MPa.
Fig. 8 Variations in the Pwave velocity under different pressures. The xaxis is the variation of porosity, and the yaxis is the Pwave velocity. The triangle symbols represent Pwave velocities at a pore pressure of 8 MPa and a CO_{2} saturation of 20%. In addition, the black colour represents the Pwave velocity under the original porosity (0.18–0.29), while the red, blue, green, grey and purple colours represent the Pwave velocities under porosities that are 1–5%, respectively, higher than the original porosity. The star symbols represent the Pwave velocities for a pore pressure of 15 MPa and a constant porosity, while the black, red, blue, and green colours represent the Pwave velocities under CO_{2} saturations of 5%, 10%, 15%, and 20%, respectively. 
The Pwave velocity varies with the porosity. However, the effects of porosity variations on the Swave velocity remain to be studied. Here, we study these effects with different pressures and CO_{2} saturations. In Figure 9, the xaxis represents the porosity, and the yaxis is the Swave velocity; the figure shows data for different pressures and saturations and porosities. Clearly, the variation in the Swave velocity is determined by the fluid pressure.
Fig. 9 Variations in the Swave velocity with the porosity. Different symbols represent different pore pressures, whereas different colours with the same symbol represent different CO_{2} saturations. 
Here, we study the changes in the Swave velocity with changes in the porosity. In Figure 10, the variation in the Swave velocity is determined mainly by the pressure, although Swave velocities at different pressures can still overlap in response to large variations in the porosity at high CO_{2} saturations (purple triangles). These results also demonstrate the significance of considering both the CO_{2} saturation and the porosity in the forward model.
Fig. 10 Variations in the Swave velocity with the porosity. The xaxis is the sample number, and the yaxis is the Swave velocity. The triangle symbols represent the Swave velocity at a pore pressure of 8 MPa and the original porosity (0.18–0.29). The black colour indicates a CO_{2} saturation of 5%, while the red, blue, green, grey, and purple colours represent CO_{2} saturations of 10%, 15%, 20%, 25%, and 30% and porosities that are 1–5%, respectively, higher than the original value. The star symbols represent a pore pressure of 15 MPa. The green, blue, red and black colours represent Swave velocities with CO_{2} saturations of 5%, 10%, 15%, and 20%, respectively, for similar porosities. 
5 Establishing a 4D seismic forward model
Based on the S and Pwave velocities obtained under different pressures, CO_{2} 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 Twolayer model
A twolayer 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.
Fig. 11 Synthetic seismograms of the twolayer model at a pore pressure of 8 MPa. The lefthand panels from left to right are the Pwave velocity, Swave velocity, density, and synthetic seismogram. The righthand panels from left to right are the Pwave velocity, Swave velocity, density, synthetic seismogram and difference profile. In the plots of the Pwave velocity, Swave velocity, and density, the curves are measured well log data, and the straight lines are block parameters. In each row, the synthetic seismogram on the far right is the difference profile between the synthetic seismograms on the left and righthand sides. a) The porosity is 0.18, while the CO_{2} saturation is 5% on the left and 20% on the right; b) The porosity is 0.24, while the CO_{2} saturation is 5% on the left and 20% on the right; c) The porosity is 0.34, while the CO_{2} saturation is 5% on the left and 20% on the right. 
The AVO gradient intercept figure (Fig. 12) corresponds to the twolayer model obtained above. As shown in Figure 12, variation in the porosity leads to variation in the AVO intercept, and variation in the CO_{2} saturation results in a change in the AVO gradient, which is consistent with the responses in the twolayer seismic model.
Fig. 12 AVO gradient intercept plot. The stars, triangles, and circles represent porosities of 0.18, 0.24, and 0.34, respectively. The red and blue colours represent CO_{2} saturations of 5% and 20%, respectively. 
To investigate the effects of porosity variations on the seismic response and AVO, synthesize seismograms for twolayer 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.
Fig. 13 Twolayer models. The figures on the left show the Pwave velocity, Swave velocity, density, and synthetic seismogram from left to right, whereas the figures on the right show the Pwave velocity, Swave velocity, density, synthetic seismogram and difference between the synthetic seismograms. In the plots of the Pwave velocity, Swave velocity and density, the curves are measured well log data, and the straight lines are block parameters. In each row, the difference between the synthetic seismograms on the left and right is shown on the far right. The conditions of the plots on the left and right of each row are as follows: a) Left: 15 MPa pore pressure, 0.2 porosity, and 5% CO_{2} saturation; right: 8 MPa pore pressure, 0.2 porosity, and 20% CO_{2} saturation; b) Left: 15 MPa pore pressure, 0.2 porosity, and 5% CO_{2} saturation; right: 8 MPa pore pressure, 0.21 porosity, and 20% CO_{2} saturation; c) Left: 15 MPa pore pressure, 0.2 porosity, and 5% CO_{2} saturation; right: 8 MPa pore pressure, 0.22 porosity, and 20% CO_{2} saturation; d) Left: 15 MPa pore pressure, 0.2 porosity, and 5% CO_{2} saturation; right: 8 MPa pore pressure, 0.23 porosity, and 20% CO_{2} saturation. 
To directly observe the effects of porosity and pressure variations on the AVO, the AVO gradient intercept is calculated based on the twolayer 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.
Fig. 14 AVO gradient intercept plot. Circles represent the conditions under a 15 MPa pore pressure, 5% CO_{2} saturation, and 0.2 porosity. Stars represent the conditions under a 15 MPa pore pressure and 20% CO_{2} saturation. Triangles represent the conditions under an 8 MPa pore pressure and 20% CO_{2} saturation. The red, blue, green, and orange colours represent porosities of 0.2, 0.21, 0.22, and 0.23, respectively. 
5.2 Well model
Based on the abovementioned twolayer models, a well model can be developed using data from Well 423 to compute synthetic seismograms and the differences between those seismograms. The well log curve from Well 423 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 Swave velocities after CO_{2} injection to fully develop the well model. Before fluid substitution, the pore pressure is 15 MPa, and the CO_{2} saturation is 0; after fluid substitution, the pore pressure is 8 MPa, and the CO_{2} 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.
Fig. 15 Synthetic seismogram of the well model. The figures on the left show the Pwave velocity, Swave velocity, density, and synthetic seismogram from left to right, whereas the figures on the right show the Pwave velocity, Swave velocity, density, synthetic seismogram and difference between the seismograms (the synthetic seismogram on the left is subtracted from that on the right). In each row, the left panel shows the original synthetic seismogram, and the right panel shows the synthetic seismogram with a consistent porosity after fluid substitution. 
6 Discussion
6.1 The effects of fluid variations on the Pwave velocity
Figure 7 shows that, at a constant pressure, the Pwave velocity decreases and then increases with increasing CO_{2} 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 Pwave velocity; therefore, with increasing CO_{2} saturation, the Pwave velocity decreases. After reaching the critical porosity, the main influence is derived from the fluid, and the variation in the Pwave 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 Pwave velocity can vary substantially at high porosity values.
Fig. 16 Properties of fluid mixtures under different pressures. The xaxis is the pore pressure. a) The yaxis is the bulk strain modulus, which varies with pressure, of the fluid mixture; b) The yaxis is the density, which varies with pressure, of the fluid mixture. 
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_{2} 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 twolayer 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 Swave 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 depthtime conversion procedure when synthesizing the seismogram.
7 Conclusion
Based on this study, we can conclude the following:

Through calculations of the P and Swave velocities of the marly layer in the Weyburn deposit under different conditions, variations in the porosity can result in variations in the P and Swave 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 Swave velocities, especially when the porosity of the deposit is relatively high.

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.

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_{2} geological storage.
References
 Avseth P., Mukerji T., Mavko G. (2006) Quantitative seismic interpretation: Applying rock physics tools to reduce interpretation risk, Cambridge University Press, Cambridge, UK. [Google Scholar]
 Berrezueta E., GonzálezMenéndeza L., Breitner D., Luquot L. (2013) Pore system changes during experimental CO_{2} injection into detritic rocks: Studies of potential storage rocks from some sedimentary basins of Spain, Int. J. Greenhouse Gas Control 17, 411–422. [CrossRef] [Google Scholar]
 Brown L.T. (2002) Integration of rock physics and reservoir simulation for the interpretation of timelapse seismic data at Weyburn Field, Saskatchewan, Master’s Thesis, Colorado School of Mines, Golden, CO. [Google Scholar]
 Digby P.J. (1981) The Effective elastic moduli of porous granular rocks, J. Appl. Mech. 48, 803–808. [Google Scholar]
 Duffaut K., Landrø M. (2007) Vp/Vs ratio versus differential stress and rock consolidation – A comparison between rock models and timelapse AVO data, Geophysics 72, 5, 81–94. [CrossRef] [Google Scholar]
 Fӧster A., Schӧner R., Fӧster H.J., Norden B., Blaschke A.W., Luckert J., Beutler G., Gaupp R., Rhede D. (2010) Reservoir characterization of a CO_{2} storage aquifer: The upper Triassic Stuttgart formation in the northeast German basin, Mar. Petrol. Geol. 27, 10, 2156–2172. [CrossRef] [Google Scholar]
 Hovorka S. (2009) Frio brine pilot: The first U.S. sequestration test, SouthWest Hydrol. 8, 26–31. [Google Scholar]
 Ivanova A., Kashubin A., Juhojuntti N., Kummerow J., Henninges J., Juhlin Ch, Lüth S., Ivandic M. (2012) Monitoring and volumetric estimation of injected CO_{2} using 4D seismic, petrophysical data, core measurements and well logging: A case study at Ketzin, Germany, Geophys. Prospect. 60, 957–973, doi: 10.1111/j.13652478.2012.01045.x. [Google Scholar]
 Jensen G.K.S. (2016) Weyburn oilfield core assessment investigating cores from pre and post CO_{2} injection: Determining the impact of CO_{2} on the reservoir, Int. J. Greenhouse Gas Control 54, Part 2, 490–498. [CrossRef] [Google Scholar]
 Karamalidis A., Hakala J.A., Griffith C., Hedges S., Lu J. (2010) Laboratory investigation of CO_{2} rockbrine interactions using natural sandstones and brine samples from the SECARB Tuscaloosa injection zone. Presented at Geological Society of America Denver Annual Meeting, 31 October – 3 November, Denver, CO, USA. [Google Scholar]
 Kharaka Y.K., Cole D.R., Hovorka S.D., Gunter W.D., Knauss K.G., Freifeld B.M. (2006) Gaswaterrock interactions in Frio Formation following CO_{2} injection: Implications for the storage of greenhouse gases in sedimentary basins, Geology 34, 577–580, doi: 10.1130/G22357.1. [Google Scholar]
 Lu J., Nicot J.P., Mickler P.J., Ribeiro L.H., Darvari R. (2016) Alteration of Bakken reservoir rock duing CO_{2}based fracturing – An autoclave reaction experiment, J. Unconv. Oil Gas Res. 14, 2016, 72–85. [CrossRef] [Google Scholar]
 Ma J., Morozov I.B. (2010) AVO modeling of pressuresaturation effects in Weyburn CO_{2} sequestration, Leading Edge 29, 2, 178–183. [CrossRef] [Google Scholar]
 Mavko G., Mukerji T., Dvorkin J.P. (1998) The rock physics handbook: 147–159, Cambridge University Press, Cambridge, UK. [Google Scholar]
 Mindlin R.D. (1949) Compliance of elastic bodies in contact, J. Appl. Mech. 16, 259–268. [Google Scholar]
 Murphy W. (1982) Effects of microstructure and pore fluids on the acoustic properties of granular sedimentary materials, PhD Thesis, Stanford University, Stanford, CA. [Google Scholar]
 Njiekak G., Schmitt D.R., Yam H., Kofman R.F. (2013) CO_{2} rock physics as part of WeyburnMidale geological storage project, Int. J. Greenhouse Gas Control 16, Suppl. 1, S118–S133. [CrossRef] [Google Scholar]
 Riding J.B., Rochelle C. (2005) The IEA Weyburn CO_{2} monitoring and storage projectfinal report of the European research team, British Geological Survey, Nottingam. [Google Scholar]
 Saul M., Lumley D. (2015) The combined effects of pressure and cementation on 4D seismic data, Geophysics 80, 2, 135–148. [CrossRef] [Google Scholar]
 Skorpa R., Todorovic J., Torsæter M. (2017) Porosity changes in mudaffected rock and cement upon reaction with CO_{2}, Energy Procedia 114, 5266–5274. [Google Scholar]
 Tambach T.J., Koenen M., Wasch L.J., Bergen F. (2015) Geochemical evaluation of CO_{2} injection and containment in a depleted gas field, Int. J. Greenhouse Gas Control 32, 2015, 61–80. [CrossRef] [Google Scholar]
 Thordsen J.J., Kharaka Y.K., Thomas R.B., Ambats G., Abedini A., Manning M.A., Lu J. (2012) Natural heterogeneity and evolving geochemistry of Lower Tuscaloosa Formation brine in response to continuing CO_{2} injection at Cranfield EOR site, Mississippi, USA, American Geophysical Union 2012 Fall Meeting, December 3–7, 2012, San Francisco, CA, United States. [Google Scholar]
 Vanorio T. (2015) Recent advances in timelapse, laboratory rock physics for the characterization and monitoring of fluidrock interactions, Geophysics 80, 2, WA49–WA59. [CrossRef] [Google Scholar]
 Wang Y., Han D.H., Ren J.L., Zhang Y., Zhao L. (2018) Microstructure effects on static and dynamic moduli for two sandstones, Meeting, SEG 2018 DPRP, 20–22 May, Beijing, China, pp. 21–24. [Google Scholar]
 White D.J. (2009) Monitoring CO_{2} storage during EOR at the WeyburnMidale Field, Leading Edge 28, 7, 838–842. [CrossRef] [Google Scholar]
 White D.J. (2013) Seismic characterization and timelapse imaging during seven years of CO_{2} flood in the Weyburn field, Saskatchewan, Canada, Int. J. Greenhouse Gas Control 16S, 78–94. [CrossRef] [Google Scholar]
All Tables
Summary of information on various porosities and their corresponding C _{ p } values.
All Figures
Fig. 1 The Phase 1A area of CO_{2} sequestration in the Weyburn field, Canada. Well 423 is the logging well, whereas Well 1213 and Well 1411 are the coring wells. 

In the text 
Fig. 2 Well log data from Well 423. From left to right are the Pwave velocity (V _{ p }), Swave velocity (V _{ s }), resistivity (Rt) log, Gamma Ray (GR) log, shale content (V _{sh}), porosity (Por; the dashed line indicates the effective porosity, while the solid line indicates the total porosity), and density. 

In the text 
Fig. 3 Linear fitting between C _{ p } and e ^{1−φ }. Black dots are Murphy’s summarization. Black line is result of linear regression. 

In the text 
Fig. 4 Comparison between the estimated and measured Swave velocities in PC1213 1410 from the marly layer. The error bar shows a 95% confidence interval for the mean of measured Swave velocity. 

In the text 
Fig. 5 Comparison between the estimated and measured Swave velocities in PC1411 1420.2 from the marly layer. The error bar shows a 95% confidence interval for the mean of measured Swave velocity. 

In the text 
Fig. 6 Variations in the Pwave velocity with different porosities under different pressure and saturation conditions. The xaxis is the porosity, and the yaxis is the Pwave velocity. The triangle symbols represent a pore pressure of 8 MPa, and the star symbols represent a pore pressure of 15 MPa. The black, green, blue, and red colours represent the Pwave velocities associated with CO_{2} saturations of 5%, 10%, 15%, and 20%, respectively. 

In the text 
Fig. 7 Variations in the Pwave velocity under different porosity, pressure and saturation conditions. The xaxis is the variation of porosity, and the yaxis is the Pwave velocity. The triangle symbols represent a pore pressure of 8 MPa, and the star symbols represent a pore pressure of 15 MPa. The black colour represents the Pwave velocity trend at a CO_{2} saturation of 5% and a porosity of 0.18–0.29, while the red, blue, green, grey and purple colours represent the Pwave velocity trends at CO_{2} saturations of 10%, 15%, 20%, 25%, and 30% and porosities that are 1–5%, respectively, higher than the data represented by the black colour. 

In the text 
Fig. 8 Variations in the Pwave velocity under different pressures. The xaxis is the variation of porosity, and the yaxis is the Pwave velocity. The triangle symbols represent Pwave velocities at a pore pressure of 8 MPa and a CO_{2} saturation of 20%. In addition, the black colour represents the Pwave velocity under the original porosity (0.18–0.29), while the red, blue, green, grey and purple colours represent the Pwave velocities under porosities that are 1–5%, respectively, higher than the original porosity. The star symbols represent the Pwave velocities for a pore pressure of 15 MPa and a constant porosity, while the black, red, blue, and green colours represent the Pwave velocities under CO_{2} saturations of 5%, 10%, 15%, and 20%, respectively. 

In the text 
Fig. 9 Variations in the Swave velocity with the porosity. Different symbols represent different pore pressures, whereas different colours with the same symbol represent different CO_{2} saturations. 

In the text 
Fig. 10 Variations in the Swave velocity with the porosity. The xaxis is the sample number, and the yaxis is the Swave velocity. The triangle symbols represent the Swave velocity at a pore pressure of 8 MPa and the original porosity (0.18–0.29). The black colour indicates a CO_{2} saturation of 5%, while the red, blue, green, grey, and purple colours represent CO_{2} saturations of 10%, 15%, 20%, 25%, and 30% and porosities that are 1–5%, respectively, higher than the original value. The star symbols represent a pore pressure of 15 MPa. The green, blue, red and black colours represent Swave velocities with CO_{2} saturations of 5%, 10%, 15%, and 20%, respectively, for similar porosities. 

In the text 
Fig. 11 Synthetic seismograms of the twolayer model at a pore pressure of 8 MPa. The lefthand panels from left to right are the Pwave velocity, Swave velocity, density, and synthetic seismogram. The righthand panels from left to right are the Pwave velocity, Swave velocity, density, synthetic seismogram and difference profile. In the plots of the Pwave velocity, Swave velocity, and density, the curves are measured well log data, and the straight lines are block parameters. In each row, the synthetic seismogram on the far right is the difference profile between the synthetic seismograms on the left and righthand sides. a) The porosity is 0.18, while the CO_{2} saturation is 5% on the left and 20% on the right; b) The porosity is 0.24, while the CO_{2} saturation is 5% on the left and 20% on the right; c) The porosity is 0.34, while the CO_{2} saturation is 5% on the left and 20% on the right. 

In the text 
Fig. 12 AVO gradient intercept plot. The stars, triangles, and circles represent porosities of 0.18, 0.24, and 0.34, respectively. The red and blue colours represent CO_{2} saturations of 5% and 20%, respectively. 

In the text 
Fig. 13 Twolayer models. The figures on the left show the Pwave velocity, Swave velocity, density, and synthetic seismogram from left to right, whereas the figures on the right show the Pwave velocity, Swave velocity, density, synthetic seismogram and difference between the synthetic seismograms. In the plots of the Pwave velocity, Swave velocity and density, the curves are measured well log data, and the straight lines are block parameters. In each row, the difference between the synthetic seismograms on the left and right is shown on the far right. The conditions of the plots on the left and right of each row are as follows: a) Left: 15 MPa pore pressure, 0.2 porosity, and 5% CO_{2} saturation; right: 8 MPa pore pressure, 0.2 porosity, and 20% CO_{2} saturation; b) Left: 15 MPa pore pressure, 0.2 porosity, and 5% CO_{2} saturation; right: 8 MPa pore pressure, 0.21 porosity, and 20% CO_{2} saturation; c) Left: 15 MPa pore pressure, 0.2 porosity, and 5% CO_{2} saturation; right: 8 MPa pore pressure, 0.22 porosity, and 20% CO_{2} saturation; d) Left: 15 MPa pore pressure, 0.2 porosity, and 5% CO_{2} saturation; right: 8 MPa pore pressure, 0.23 porosity, and 20% CO_{2} saturation. 

In the text 
Fig. 14 AVO gradient intercept plot. Circles represent the conditions under a 15 MPa pore pressure, 5% CO_{2} saturation, and 0.2 porosity. Stars represent the conditions under a 15 MPa pore pressure and 20% CO_{2} saturation. Triangles represent the conditions under an 8 MPa pore pressure and 20% CO_{2} saturation. The red, blue, green, and orange colours represent porosities of 0.2, 0.21, 0.22, and 0.23, respectively. 

In the text 
Fig. 15 Synthetic seismogram of the well model. The figures on the left show the Pwave velocity, Swave velocity, density, and synthetic seismogram from left to right, whereas the figures on the right show the Pwave velocity, Swave velocity, density, synthetic seismogram and difference between the seismograms (the synthetic seismogram on the left is subtracted from that on the right). In each row, the left panel shows the original synthetic seismogram, and the right panel shows the synthetic seismogram with a consistent porosity after fluid substitution. 

In the text 
Fig. 16 Properties of fluid mixtures under different pressures. The xaxis is the pore pressure. a) The yaxis is the bulk strain modulus, which varies with pressure, of the fluid mixture; b) The yaxis is the density, which varies with pressure, of the fluid mixture. 

In the text 