The AVO Effect of Formation Pressure on TimeLapse Seismic Monitoring in Marine Carbon Dioxide Storage
https://doi.org/10.1007/s1180402400439w

Abstract
The phase change of CO_{2} has a significant bearing on the siting, injection, and monitoring of storage. The phase state of CO_{2} is closely related to pressure. In the process of seismic exploration, the information of formation pressure can be response in the seismic data. Therefore, it is possible to monitor the formation pressure using timelapse seismic method. Apart from formation pressure, the information of porosity and CO_{2} saturation can be reflected in the seismic data. Here, based on the actual situation of the work area, a rockphysical model is proposed to address the feasibility of timelapse seismic monitoring during CO_{2} storage in the anisotropic formation. The model takes into account the formation pressure, variety minerals composition, fracture, fluid inhomogeneous distribution, and anisotropy caused by horizontal layering of rock layers (or oriented alignment of minerals). From the proposed rockphysical model and the welllogging, cores and geological data at the target layer, the variation of Pwave and Swave velocity with formation pressure after CO_{2} injection is calculated. And so are the effects of porosity and CO_{2} saturation. Finally, from anisotropic exact reflection coefficient equation, the reflection coefficients under different formation pressures are calculated. It is proved that the reflection coefficient varies with pressure. Compared with CO_{2} saturation, the pressure has a greater effect on the reflection coefficient. Through the convolution model, the seismic record is calculated. The seismic record shows the difference with different formation pressure. At present, in the marine CO_{2} sequestration monitoring domain, there is no study involving the effect of formation pressure changes on seismic records in seafloor anisotropic formation. This study can provide a basis for the inversion of reservoir parameters in anisotropic seafloor CO_{2} reservoirs.Article Highlights● A seafloor anisotropic rockphysical model that considering formation pressure is proposed.● The anisotropic exact reflection coefficient equation is used in seismic forward modeling.● The influence of formation pressure on velocity, reflection coefficient, and seismic record is clarified.● In addition to formation pressure, the saturation and porosity of CO_{2} also have different effects on seismic records. 
1 Introduction
Subsurface rockophysical modeling and timelapse seismic monitoring are two key technologies in CO_{2} storage. The former can establish the bridge between elastic and reservoir parameters by simulating the actual situation of the seabed formation; the latter is to obtain the changes of subsurface velocity, density, impedance and other parameters through inversion, so as to evaluate the location and state of CO_{2}. Most CO_{2} sequestration projects worldwide involve rockophysics and timelapse seismic monitoring. Among them, several famous projects include the WeyburnMidale project in Canada (IEA GHG WeyburnMidale CO_{2} Monitoring and Storage project) (White, 2009, 2013; Smith et al., 2018), the In Salah project of BP in Algeria, the Aquistore project of Saskatchewan Power Plant and Petroleum Technology Research Centre in Canada, and the Ketzin Project of the European Union in Germany. In terms of offshore CO_{2} storage, the famous projects that have been in operation are the Sleipner and Snovit project in Norway, and the Tomakomai project in Japan. In Europe, there are several offshore carbon storage projects centered in the North Sea. In the future, the Northern Lights in Norway and Greensand in Denmark, are expected to be active in 2024 and 2025, respectively (Patil et al., 2023).
Although rockphysical modeling and timelapse seismic monitoring have been successful operated in many domestic and international CO_{2} sequestration projects (Alfi and Hosseini, 2016; Wang et al., 2017), the case has not been applied to offshore carbon sequestration in China. The temperature and pressure conditions in the shallow seabed formation, the seismic data acquisition is different from land condition. And there are also many differences in implementation and expense. At present, there is a lack of expense research and understanding of whether CO_{2} offshore sequestration can be carried out with better timelapse seismic monitoring effects and the physical significance of the differences in timelapse seismic monitoring. After the current smallscale experimental CO_{2} injection in the South China Sea in China, the simulation study of its geology environment and the forward analysis of timelapse seismic monitoring, especially the effect of formation pressure on seismic data, need to be further studied.
The establishment of a timelapse seismic forward model is the basis of seismic inversion. The main factor for accurate calculation of seismic forward modeling is the elastic parameters during CO_{2} injection. Especially the Pwave and Swave velocity and density. Because they are the most prominent parameters affecting the timelapse seismic forward modeling. During the CO_{2} injection process, the reservoir pressure as well as fluid saturation changed, which affected the Pwave and Swave velocity and density. In the offshore, the cost of drilling and logging is very high. Furthermore, the reservoir information acquired by logging methods is very limited. Therefore, the detailed rockphysical model is necessary to simulate the reservoir. Then, with the help of anisotropic exact reflection coefficient equation, the more accurate forward can be obtained (Ma and Morozov, 2010).
In the situation that no events of fracturing, collapse, and subsidence before and after CO_{2} injection, the difference in amplitude and travel time between the two seismic monitoring mainly comes from the effect of fluid saturation and formation pressure (Li et al., 2017). This condition is the result of reservoir elastic parameters varied by saturation and pressure. Understanding the changing trend of elastic parameters with formation pressure and CO_{2} saturation is important to timelapse seismic forward modeling. Another way to obtaining this relationship is cores analyzing. It can get a statistical relationship. However, such statistical results are limited by the quantity and location of cores. And this method lacks the support of rockphysical theories (Wang et al., 2019, Wang et al., 2022). This lead to a poor generalization. Therefore, it is very necessary to establish a more reasonable and comprehensive research method based on rockphysical theory (Wang et al., 2020).
Specifically, in conjunction with the actual seafloor formation geological situation, a rockphysical model that can take into account varieties of influencing factors (mineral composition, pore aspect ratio, fracture, fluid distribution, formation pressure, and anisotropy) is proposed. Through the rockphysical model, the change of Pwave and Swave velocity during CO_{2} injection can be obtained. Then the reflection coefficient is calculated by the anisotropic exact reflection coefficient equation. With the reflection coefficients, together with the commonly used Ricker wavelet, the seismic record is calculated by the convolution model. In this way, the relationship between formation pressure (and CO_{2} saturation) and seismic record is established. This contributes to analyzing the influence of formation pressure and CO_{2} saturation on the timelapse amplitude versus offset (AVO) interpretation during the CO_{2} injection process. The conclusion provides a theoretical basis for the timelapse seismic interpretation and monitoring.
2 Rockphysics modeling
2.1 Equivalent elastic modulus of VTI background medium
For rocks with complex mineral compositions, the equivalent medium theory can be used to calculate the modulus of mixed minerals. Compared with the frequently used VoigtReuss boundary and Hill's average, the upper and lower bounds of the Hashin Shtrikman boundary (Hashin and Shtrikman, 1962, 1963) are the narrowest allowable upper and lower bounds (Mavko et al., 2020). Therefore, the arithmetic average value of HashinShtrikman boundary is used as the equivalent modulus of mixed minerals. The original HashinShtrikman boundary is suitable for two components. Berryman (1995) deduced a more general form of HashinShtrikman boundary for calculating more than two components:
$$ \begin{aligned} & K^{\mathrm{HS}+}=\Lambda\left(\frac{4 \mu_{\max }}{3}\right), K^{\mathrm{HS}}=\Lambda\left(\frac{4 \mu_{\min }}{3}\right) \\ & \mu^{\mathrm{HS}+}=\Gamma\left(\zeta\left(K_{\max }, \mu_{\max }\right)\right), \mu^{\mathrm{HS}}=\Gamma\left(\zeta\left(K_{\min }, \mu_{\min }\right)\right) \end{aligned} $$ (1) where,
$$ \begin{aligned} \Lambda(z) & =\left\langle\frac{1}{K(r)+z}\right\rangle^{1}z, \Gamma(z)=\left\langle\frac{1}{\mu(r)+z}\right\rangle^{1} \\ & z, \zeta(K, \mu)=\frac{\mu}{6}\left(\frac{9 K+8 \mu}{K+2 \mu}\right) \end{aligned} $$ (2) K^{HS+} and K^{HS}^{−} are the upper and lower bounds of bulk modulus, μ^{HS+} and μ^{HS}^{−} are the upper and lower bounds of shear modulus, K(r) and μ(r) are bulk modulus and shear modulus of each component. The braces indicate the average of materials, that is, the weighted average of each component according to its volume content.
The HertzMindlin contact theory is to equate the rock as an accumulation of some samesized grains, by calculating the positive and tangential stiffnesses between two grains under certain external forces, and then obtaining the equivalent modulus of an arbitrary equalsize accumulation of grains (Mindlin, 1949). On the basis of the HertzMindlin contact theory, many scholars have made different aspect improvements to this theory (Walton, 1987; Digby, 1981; Jenkins et al., 2005; Norris and Johnson, 1997; Dvorkin and Nur, 1996). Among them, Dvorkin improved the model by taking into account the sorting of the grains inside the rock. In this way, though considering the relaxation between grains and the irregularity of the contact surface among grains, it is possible to simulate the characteristics of the actual seafloor formation. Especially the weak consolidation and variation of grains sizes, etc. The expression of Dvorkin's improvement can be summarized in the following form.
$$ K_{\mathrm{eff}}=\left(\frac{(1\phi)^2 \mu^2 n^2 P}{18 \pi^2(1v)^2} \cdot \frac{\bar{R}}{R}\right)^{\frac{1}{3}} $$ (3) $$ \begin{aligned} \mu_{\mathrm{eff}}= & \frac{1}{10}\left(\frac{12(1\phi)^2 \mu^2 n^2 P}{\pi^2(1v)^2} \cdot \frac{\bar{R}}{R}\right)^{\frac{1}{3}} \\ & +\frac{3 f_t}{10}\left(\frac{12(1\phi)^2 \mu^2 n^2 P(1v)}{\pi^2(2v)^3} \cdot \frac{\bar{R}}{R}\right)^{\frac{1}{3}} \end{aligned} $$ (4) where, ν is the Poisson's ratio of mineral, R is the curvature radius of the contact surface, f is nonsmooth contact percentage, R is the effective contact radius, P is the effective pressure at the point of grains contact. The deduction in HertzMindlin theory involves the HashinShtrikman boundary model(Hashin and Shtrikman, 1962, 1963). The equation can be used to calculate the bulk and shear modulus at a given depth by assuming mineral composition, porosity, and overburden stress. As shown in equations (3) and (4), unless the granular aggregates are composed of the same spherical grains (R/R=1), the calibration parameter are dependent on the coordination number n and ratio R/R. R/R can be greater than or less than 1. The ratio takes into account the smoothness and sorting of grains, and can be regarded as a parameter that characterizes the average effective contact radius and general grains size distribution.
After obtaining the modulus of various mineral mixtures, the stiffness coefficient of the VTI background medium can be calculated. The elastic parameters of anisotropic thin layered strata composed of multiple isotropic layers are calculated by Backus average (Backus, 1962). Although the assumption of Backus average is long wavelength, for horizontally layered anisotropic medium, the usage on the logging scale is effective (Hsu et al., 1988; Imhof, 2003; Liner, 2006). The VTI equivalent medium by Backus average can be written as the stiffness matrix,
$$ \boldsymbol{c}_{\mathrm{VTI}}^B=\left[\begin{array}{cccccc} c_{11}^B & c_{12}^B & c_{13}^B & 0 & 0 & 0 \\ c_{12}^B & c_{11}^B & c_{13}^B & 0 & 0 & 0 \\ c_{13}^B & c_{13}^B & c_{33}^B & 0 & 0 & 0 \\ 0 & 0 & 0 & c_{44}^B & 0 & 0 \\ 0 & 0 & 0 & 0 & c_{44}^B & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{c_{11}^Bc_{12}^B}{2} \end{array}\right] $$ (5) where, the relationship between stiffness coefficient and elastic parameters of each layer is as follows:
$$ \begin{aligned} & c_{11}^B=\left\langle\frac{4 \mu(\lambda+\mu)}{\lambda+2 \mu}\right\rangle+\left\langle\frac{1}{\lambda+2 \mu}\right\rangle^{1}\left\langle\frac{\lambda}{\lambda+2 \mu}\right\rangle^2 \\ & c_{12}^B=\left\langle\frac{2 \mu \lambda}{\lambda+2 \mu}\right\rangle+\left\langle\frac{1}{\lambda+2 \mu}\right\rangle^{1}\left\langle\frac{\lambda}{\lambda+2 \mu}\right\rangle^2 \\ & c_{13}^B=\left\langle\frac{1}{\lambda+2 \mu}\right\rangle^{1}\left\langle\frac{\lambda}{\lambda+2 \mu}\right\rangle \\ & c_{33}^B=\left\langle\frac{1}{\lambda+2 \mu}\right\rangle^{1} \\ & c_{55}^B=\left\langle\frac{1}{\mu}\right\rangle^{1} \\ & c_{66}^B=\langle\mu\rangle \end{aligned} $$ (6) where superscript B represents the calculation result of Backus average, λ and μ are Lamé coefficient.
2.2 Equivalent elastic modulus of dry rock skeleton with pores
In the seafloor, the porosity is relative high. The effective medium theory can be used for simulation this condition. Hornby et al. (1994) studied the effective elastic properties by using the effective medium model. The model combines anisotropic selfconsistent approximations model (SCA) with anisotropic differential effective medium model (DEM). Many scholars (Yuan, 2007; Qian, 2017; Zhao, 2017) have obtained two perceptions after forward simulation calculations. One is that the anisotropic SCA can ensure that the matrix and fluid of rock are connected in twophases only when the porosity is within 40%‒60%. Therefore, for the solid – liquid twophase mixture, to simulate the actual twophase connection rock structure, the anisotropic SCA model is used to calculate the equivalent stiffness tensor when the fluid content is 50%, and then the anisotropic DEM model is used to adjust to the specific porosity to obtain the calculation results of anisotropic SCA+DEM. The other is a disadvantage of DEM. The DEM formulation does not treat each constituent symmetrically. There is a preferred matrix or host material, and the effective moduli depend on the construction path taken to reach the final composite.
For multiple inclusion shapes or multiple constituents, the effective moduli depend not only on the final volume fractions of the constituents but also on the order in which the incremental additions are made. The anisotropic SCA+ DEM model can be written as.
$$ \frac{\mathrm{d}}{\mathrm{~d} v}(c(v))=\frac{1}{1v}\left(c^ic(v)\right) \cdot\left[I+\hat{G}\left(c^ic(v)\right)\right]^{1} $$ (7) where, c represents the calculated equivalent medium stiffness tensor; v is the volume content of the inclusion; c^{i} is the stiffness tensor of the inclusion; The tensor $\stackrel{\wedge}{G}$ is Eshelby tensor (Eshelby, 1957), which controls the shape of pores, and its expression contains the aspect ratio. The derivation and calculation process of $\stackrel{\wedge}{G}$ are given in the form of line integral (Kinoshita and Mura, 1971; Mura, 1991).
2.3 Equivalent elastic modulus of dry rock skeleton with fractures
Overpressure fluids and gases can cause the formation of fractures in muddy sediments (Sassen et al. 2001). Figure 1 shows a fracture filled natural gas hydrate model, Using the linear slip theory, the VTI background medium is embedded with horizontal direction parallel fractures (Hsu et al. 1988; Sondergeld and Rai, 2011). With the help of this step, the rockphysical model considers both the anisotropy caused by the directional arrangement of minerals and fractures. This is more relevant to the reality of the underground situation. The assumptions of the model are consistent with Schoenberg and Helbig's model, ignoring the influence relationship between horizontal and vertical symmetry. The elastic coefficient matrix form of the above theory can be written as follows:
$$ \boldsymbol{C}_{\mathrm{VTI}}=\left[\begin{array}{cccccc} c_{11 b}(1\Delta N) & c_{12 b}(1\Delta N) & c_{13 b}(1\Delta N) & 0 & 0 & 0 \\ c_{12 b}(1\Delta N) & c_{11 b}\Delta N \frac{c_{12 b}}{c_{11 b}} & c_{13 b}\left(1\Delta N \frac{c_{12 b}}{c_{11 b}}\right) & 0 & 0 & 0 \\ c_{13 b}(1\Delta N) & c_{13 b}\left(1\Delta N \frac{c_{12 b}}{c_{11 b}}\right) & c_{33 b}\Delta N \frac{c_{13 b}}{c_{11 b}} & 0 & 0 & 0 \\ 0 & 0 & 0 & c_{44 b} & 0 & 0 \\ 0 & 0 & 0 & 0 & c_{44 b}(1\Delta T) & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{c_{11 b}c_{12 b}}{2}(1\Delta T) \end{array}\right] $$ (8) where subscript b represents the stiffness coefficient of VTI background matrix. Combining Schoenberg and Helbig's theory with Hudson's theory (Hudson, 1980, 1981), the expressions of normal weakness ΔN and tangential weakness ΔT can be obtained based on linear anisotropy approximation.
$$ \begin{aligned} & \Delta N=\frac{4 e_d}{3 g(1g)\left[1+\frac{1}{\pi(1g)} \frac{4 k^{\prime}+3 \mu^{\prime}}{4 \mu} \frac{a_l}{c_w}\right]} \\ & \Delta T=\frac{16 e_d}{3 g(32 g)\left[1+\frac{4}{\pi(32 g)} \frac{\mu^{\prime}}{\mu} \frac{a_l}{c_w}\right]} \end{aligned} $$ (9) where, g = V_{S}^{2}/V_{P}^{2}, ε^{V} and δ^{V} are two of anisotropic parameters proposed by Rüger (1997) on the analogy with Thomsen VTI anisotropic coefficients (Thomsen, 1986):
$$ \begin{gathered} \varepsilon^V=\frac{c_{11}c_{33}}{2 c_{33}}, \gamma^V=\frac{c_{66}c_{44}}{2 c_{44}} \\ \delta^V=\frac{\left(c_{13}+c_{44}\right)^2\left(c_{33}c_{44}\right)^2}{2 c_{33}\left(c_{33}c_{44}\right)} \end{gathered} $$ (10) 2.4 Calculation of mixed fluid modulus and fluid substitution
For saturated rocks with two or more fluids, the equivalent bulk modulus of the mixed fluids is required for fluid replacement. The earliest theory on fluid mixing modulus was proposed by Wood and Hill as the upper and lower limits (Wood, 1955; Hill, 1963). For the case of gaswater saturation, Brie et al. (1995) proposed an equivalent bulk modulus empirical equation. The above is calculating the mixed fluid modulus from the perspective of fluid distribution. From the aspect of rock microstructural description, harmonic and arithmetic averaging at the endpoints and capillary pressure in pores, the mixed fluid modulus is described as (Papageorgiou et al., 2016):
$$ K_q=\frac{S_w K_w+\alpha\left(1S_w\right) K_g}{S_w+\alpha\left(1S_w\right)}, 1 \leqslant \alpha \leqslant \frac{K_w}{K_g} $$ (11) This equation has a solid theoretical foundation. The theory exactly produces both a serial and a parallel law. For the limiting values of the capillary pressure coefficient parameter α, equation (11) reduces analytically to harmonic (α = K_{w}/K_{g}) or arithmetic (α = 1) averaging at the endpoints.
The BrownKorringa theory (Brown and Korringa, 1975) describes the relationship between the effective elastic tensor of anisotropic dry rock and fluidsaturated rock. It is commonly regarded as the improved anisotropic Gassmann model. Its equation can be expressed as equation (12).
$$ S_{\mathrm{ijkl}}^{\mathrm{sat}}=S_{\mathrm{ijkl}}^{\mathrm{dry}}\frac{\left(S_{\mathrm{ij} \alpha a}^{\mathrm{dry}}S_{\mathrm{ij} \alpha a}^0\right)\left(S_{\mathrm{kl} \alpha a}^{\mathrm{dry}}S_{\mathrm{kl} \alpha a}^0\right)}{\left(S_{a \alpha \beta \beta}^{\mathrm{dry}}S_{\alpha a \beta \beta}^0\right)+\left(\rho_y\rho_0\right) \phi} $$ (12) where S_{ijkl}^{dry}, S_{ijkl}^{sat}, and S_{ij}_{αα}^{0} represent the effective elastic compliance of dry rock, fluid saturated rock and solid minerals respectively. β_{fl} = 1/K_{fl} and β_{0} = S_{a}_{αββ}^{0} = 1/K_{0} represent the fluid and mineral compressibility respectively. ϕ is rock porosity.
At last but not least, the vertical direction Pwave velocity and Swave can be calculated by:
$$ V_P=\sqrt{\frac{c_{33}}{\rho}}, V_S=\sqrt{\frac{c_{44}}{\rho}} $$ (13) The flow chart of the rockphysical modeling process is shown in Figure 1 and Figure 2.
3 The effect of formation pressure, porosity and CO_{2} saturation on velocity
From the rockphysical model proposed above, the effect of formation pressure, porosity and CO_{2} saturation on the velocity can be obtained by inputting information and data from welling, logging, and cores in an area in the South China Sea.
The actual data comes from the Pearl River Estuary Basin. The Pearl River Estuary Basin is a Cenozoic continental margin extensional basin formed on the Tertiary basement, and is also a major oil and gas producing area in South China Sea. The specific location is on the west side of the depression. This region is rich in material sources and has a complex and diverse composition of minerals. Multiple wells have produced oil and gas flows since 2016, revealing the enormous exploration potential of the region. Due to multiple factors such as sedimentary environment, diagenesis, secondary transformation of oil and gas migration, and engineering drilling, this area currently has the following characteristics: shallow burial (1 000‒1 400 m), thin interbed layers (1‒4 m), fine layer grains (argillaceous siltstone, thin interbedded siltstone), and locally developed horizontal fractures. These features are consistent with the characteristics of proposed rockphysical model.
The parameters of target layer are shown in the table, which are the input for rockphysical forward modeling.
Minerals Bulk modulus (GPa) Shear modulus (GPa) Density (g/cm^{3}) Quartz 36 45 2.65 Felspar 38 17 2.68 Calcite 71 32 2.70 Clay 19 9 2.43 Brine 2.49 ‒ 1.06 CO_{2} 2.32 ‒ 0.96 The first is the effect of formation pressure on the velocity. Because the injected CO_{2} pressure is greater than the formation pressure, the CO_{2} injection will definitely affect the formation pressure, which in turn changing the Pwave and Swave velocity. In Figure 3, the formation pressure is positively correlated with the velocity. Analyzed from the geological point of view, the increase of formation pressure mainly increases the degree of compaction of rock skeleton and increases the modulus, which in turn makes the velocity increase. This can also lead to pore pressure increasing. Especially existing the overlying gap layer, the fluid space is relatively limited, and the increase in pore pressure may be more pronounced. Moreover, the elastic parameters in supercritical state CO_{2} vary more with pressure than in seawater.
The effect of porosity on velocity can also be concluded from the rockphysical model. In Figure 4, it can be seen that the porosity is negatively related to the velocity. One of the main reasons for this is that the formation become harder, denser, modulus getting bigger, which can cause velocity increasing. In addition, after CO_{2} injection, under suitable temperature and pressure conditions, CO_{2} hydrate may be generated within the formation, resulting in porosity decreasing, which will also lead to the velocity increasing.
The last important influence on CO_{2} sequestration is the saturation of CO_{2}. In Figure 5, as CO_{2} injected, the saturation of CO_{2} becomes higher, the saturation of seawater becomes lower, the Pwave velocity decreases, and the Swave velocity remains essentially unchanged. The Swave velocity is mainly affected by the solid medium, so the change of the type and amount of fluid in the formation pore space has almost no effect on the Swave. Relatively, the impact of CO_{2} saturation on velocity is not as significant as the impact of formation pressure on velocity. It is a bit difficult to detect saturation using seismic methods.
Through the forward analysis of the rockphysical model above, the formation pressure has the greatest impact on velocity, followed by porosity and CO_{2} saturation. Thus, the formation pressure has a significant impact on the elastic properties of relatively loose seabed formation.
4 The effect of formation pressure, porosity and CO_{2} saturation on reflection coefficient
Formation pressure, porosity and CO_{2} saturation have an effect on the velocity, so these two parameters certainly have an effect on the reflection coefficient. Based on the anisotropic exact reflection coefficient equation (Wu et al., 2022), it is possible to simulate the seafloor formation and obtain the variation of the reflection coefficient. This reflection coefficient equation corresponds to the rockphysical model proposed above, applicable to anisotropic formation. The equation is not a linear approximate equation and have high computational accuracy.
As shown in Figure 6, varieties of formation pressure, porosity, and CO_{2} saturation cause different reflection coefficient. The changing of reflection coefficient caused by formation pressure is the most pronounced, especially for big incidence angle. Porosity also has an effect on the reflection coefficient. The bigger the porosity, the smaller the reflection coefficient. This is in accordance with general geologic knowledge. For the effect of CO_{2} saturation on reflection coefficient, as mentioned above in the manuscript, changes in the CO_{2} content in the pore and fracture space do not affect the Swave velocity and have little effect on the Pwave velocity. There is also little difference in density between supercritical CO_{2} and seawater. This gives rise to the fact that changes in the CO_{2} saturation have almost no effect on the reflection coefficient. That is, the reflection coefficient basically only varies with incidence angle (known as the AVO effect), as shown in Figure 6(c).
5 The effect of formation pressure, porosity and CO_{2} saturation on seismic record
AVO method is used to study the seismic reflection amplitude variation characteristics between the shot point and the receiver the shot offset, and then detect the underlying media information at the reflection interface. AVO is currently the most important and commonly used method for reservoir exploration.
The VTI media exact reflection coefficient equation has no assumptions of small incident angle and weak impedance differences. So it is more applicable and computational accuracy (Wu et al. 2022). The specific form of the equation can be found in Appendix. Then using convolutional model, VTI media seismic records are obtained from the calculated reflection coefficients. In the calculation, we use the commonly used Ricker wavelet with a main frequency of 30Hz. The incidence angle is 6°‒ 26°, with an interval of 3°.
During the CO_{2} injection stage, most logging methods cannot be used due to the influence of casing in the well. Moreover, the CO_{2} saturation and formation pressure changed simultaneously. Here, reservoir pressure data is obtained through pressure sensors on the wellbore at the perforation location. The saturation of CO_{2} in the reservoir can be roughly calculated based on parameters such as the cumulative injection of CO_{2}, reservoir capacity, and porosity. In summary, through calculation, during the process of CO_{2} saturation ranging from 20% to 60%, the formation pressure increased from 1.2 MPa to 1.6 MPa. The seismic records performed before and after this process are shown in Figures 7(a) ‒ (b). And the difference between the two seismic records is shown in Figure 7(c). Figure 7(c) illustrate the difference in seismic records mainly occur in far offset and large wave impedance differences areas. This is because formation pressure has a greater impact on velocity changes at large offset distance (incidence angle). And big wave impedance differences in this area could increase this impact. The target layer is located near 0.51 s twoway travel time. In Figure 7(b), the AVO characteristic presents a strong negative polarity amplitude with zero offset and have an increasing trend. This layer has a lower wave impedance than the overlying medium. This characteristic is very similar to the class Ⅲ AVO characteristic. Class Ⅲ AVO characteristic manifested as a poorly compacted or unconsolidated low impedance gas bearing sandstones reservoir characteristic (Ismail et al., 2020). This also indicates that the seismic response of supercritical CO_{2} is similar to natural gas. Fortunately, the class Ⅲ AVO characteristic can easily detected in shallow seabed formations by conventional seismic techniques.
6 Conclusions
The construction difficulty, period and cost of offshore seismic surveys are much larger than onshore condition. Currently, there are many successful cases of onshore carbon sequestration projects in China. However, the domestic offshore carbon sequestration project has just been carried out, so applying offshore seismic survey to image the CO_{2} sequestration reservoir requires lots of work. Here, a rockphysical model that takes into account the various minerals composition of rock matrix and various fluids composition, fracture, fluid distribution, and anisotropy is proposed for the seabed carbon sequestration area. The submarine formation is simulated in sufficient detail and determine the relationship between the elastic and reservoir parameters. After that, the reflection coefficient is calculated from the anisotropic exact reflection coefficient equation, and the seismic record is obtained from the convolution model. After the above steps, it is possible to forward modeling the effect of reservoir parameters (especially pressure, the main parameter affecting the phase state of CO_{2} in shallow submarine formation) on the seismic record. The three forward modeling (physical parameters on velocity, reflection coefficients and seismic records) lead to the following insights:
1) Rockphysical modeling for the target layer is the foundation of timelapse seismic monitoring. The establishment of a rockphysical model can simulate submarine geological characteristics in detail and introduce reservoir parameters in seismic domain.
2) After filling CO_{2} into the pore and fracture, the formation pressure has a greater effect on the velocity, reflection coefficient, and seismic record. Specifically, it brings about the velocity and reflection coefficients increasing.
3) After injecting CO_{2} into the formation, the "bright spot" that similar to the Class Ⅲ AVO characteristic in unconsolidated formation may be shown in the seismic record. This indicates that it is possible to carry out timelapse seismic monitoring at different CO_{2} marine geological storage stages.
Appendix An exact reflection coefficient equation for VTI media
The Thomsen anisotropy parameters (Thomsen, 1986) mainly describe the relationship between two vertical velocities, three dimensionless anisotropy parameters and five stiffness coefficients based on the weak anisotropy assumption. The relations are as follows:
$$ \begin{aligned} & V_P=\sqrt{\frac{c_{33}}{\rho}}, V_S=\sqrt{\frac{c_{55}}{\rho}}, \varepsilon=\frac{c_{11}c_{33}}{2 c_{33}} \\ & \gamma=\frac{c_{66}c_{44}}{2 c_{55}}, \delta=\frac{\left(c_{13}+c_{55}\right)^2\left(c_{33}c_{55}\right)^2}{2 c_{33}\left(c_{33}c_{55}\right)} \end{aligned} $$ (A1) where V_{P} and V_{S} denote the P and S wave velocity in the vertical direction (symmetry axis direction) respectively, ρ denotes the density, and ε, δ, and γ are the Thomsen anisotropy parameters. Since the SHwave is completely decoupled from the P and SV waves in VTI media, the reflection and transmission coefficients of the P and SV waves are considered.
Based on the weak anisotropy assumption, the Pwave phase velocity and stiffness coefficient c_{13} in the VTI media can be linearly expressed by anisotropy parameters ε and δ as (Thomsen, 1986).
$$ V_P(\theta) \approx V_{P 0}\left(1+\delta \sin ^2 \theta \cos ^2 \theta+\varepsilon \sin ^4 \theta\right) $$ (A2) $$ c_{13} \approx c_{33}(1+\delta)c_{55} $$ (A3) According to Equation (A2), a relatively simple form of horizontal slowness (ray parameter) can be obtained.
$$ p=\frac{\sin \theta}{V_P\left(1+\delta \sin ^2 \theta \cos ^2 \theta+\varepsilon \sin ^4 \theta\right)} $$ (A4) Substituting Equation (A1) and (A3) into the VTI media exact reflection coefficient equation (Graebner, 1992; Schoenberg and Protázio, 1992; Rüger, 1996). Then the following matrix can be obtained:
$$ \left[\begin{array}{llll} m_{11} & m_{12} & m_{13} & m_{14} \\ m_{21} & m_{22} & m_{23} & m_{24} \\ m_{31} & m_{32} & m_{33} & m_{34} \\ m_{41} & m_{42} & m_{43} & m_{44} \end{array}\right]\left[\begin{array}{l} R_{\mathrm{PP}}^{\mathrm{VTI}} \\ R_{\mathrm{PS}}^{\mathrm{VTI}} \\ T_{\mathrm{PP}}^{\mathrm{VTI}} \\ T_{\mathrm{PS}}^{\mathrm{VTI}} \end{array}\right]=\left[\begin{array}{c} m_{11} \\ m_{21} \\ m_{31} \\ m_{41} \end{array}\right] $$ (A5) Like exact Zoeppritz equation (Zoeppritz, 1919; Aki and Richards, 1980), equation (A5) can be expressed in the following equation:
$$ M \cdot S=N $$ (A6) In Equation (A5), the coefficients of the matrix are given by,
(A7) The superscripts ^{(1)} and ^{(2)} indicate the upper and lower layer respectively,
where,
$$ \begin{aligned} & l_\alpha=\sqrt{\frac{V_{P 0}^2 q_\alpha^2+V_{S 0}^2 p^21}{(2 \varepsilon+1) V_{P 0}^2 p^2+V_{S 0}^2 q_\alpha^2+V_{P 0}^2 q_\alpha^2+V_{S 0}^2 p^22}} \\ & m_\alpha=\sqrt{\frac{V_{S 0}^2 q_\alpha^2+(2 \varepsilon+1) V_{P 0}^2 p^21}{(2 \varepsilon+1) V_{P 0}^2 p^2+V_{S 0}^2 q_\alpha^2+V_{P 0}^2 q_\alpha^2+V_{S 0}^2 p^22}} \\ & l_\beta=\sqrt{\frac{V_{S 0}^2 q_\beta^2+(2 \varepsilon+1) V_{P 0}^2 p^21}{(2 \varepsilon+1) V_{P 0}^2 p^2+V_{S 0}^2 q_\beta^2+V_{P 0}^2 q_\beta^2+V_{S 0}^2 p^22}} \\ & m_\beta=\sqrt{\frac{V_{P 0}^2 q_\beta^2+V_{S 0}^2 p^21}{(2 \varepsilon+1) V_{P 0}^2 p^2+V_{S 0}^2 q_\beta^2+V_{P 0}^2 q_\beta^2+V_{S 0}^2 p^22}} \end{aligned} $$ (A8) In the above equations for the VTI reflection and transmission coefficients, the vertical slowness of the different wave modes can be expressed as (Rüger, 1996):
$$ q_\alpha=\frac{\sqrt{K_1\sqrt{K_1^24 K_2 K_3}}}{\sqrt{2}}, q_\beta=\frac{\sqrt{K_1+\sqrt{K_1^24 K_2 K_3}}}{\sqrt{2}} $$ (A9) where,
$$ \begin{aligned} & K_1=\frac{\rho}{c_{33}}+\frac{\rho}{c_{55}}\left(\frac{c_{11}}{c_{55}}+\frac{c_{55}}{c_{33}}\frac{\left(c_{13}+c_{55}\right)^2}{c_{33} c_{55}}\right) p^2 \\ & K_2=\frac{c_{11}}{c_{33}} p^2\frac{\rho}{c_{33}}, K_3=p^2\frac{\rho}{c_{55}} \end{aligned} $$ (A10) and q_{α} is the vertical slowness of the Pwave, q_{β} is the vertical slowness of the SVwave, c_{ij} is the stiffness tensor, and θ is the phase angle. The subscript α and β denote the Pwave mode and SVwave mode respectively.
Substituting Equations (A1) and (A2) into Equation (A10), replaced the independent stiffness coefficients by the Thomsen anisotropy parameters.
Equation (A5) is an implicit equation of the VTI media reflection and transmission coefficients. Equations (A4), (A7), (A8), and (A9) are the intermediate variables of Equation (A5).
Competing interest The authors have no competing interests to declare that are relevant to the content of this article. 
Table 1 The parameters of target layer
Minerals Bulk modulus (GPa) Shear modulus (GPa) Density (g/cm^{3}) Quartz 36 45 2.65 Felspar 38 17 2.68 Calcite 71 32 2.70 Clay 19 9 2.43 Brine 2.49 ‒ 1.06 CO_{2} 2.32 ‒ 0.96 
Aki K, Richards PG (1980) Quantitative seismology: Theory and methods. San Francisco: W. H. Freeman, 1–932 Alfi M, Hosseini SA (2016) Integration of reservoir simulation, history matching, and 4D seismic for CO_{2}EOR and storage at Cranfield. Fuel 175: 116–128. https://doi.org/10.1016/j.fuel.2016.02.032 Backus GE (1962) Longwave elastic anisotropy produced by horizontal layering. Journal of Geophysical Research 149(67): 4427–4440. https://doi.org/10.1029/JZ067i011p04427 Berryman JG (1995) Mixture theories for rock properties. American Geophysical Union 3: 205–228. https://doi.org/10.1029/RF003p0205 Brie A, Pampuri F, Marsala AF, Meazza O (1995) Shear sonic interpretation in gasbearing sands. SPE Annual Technical Conference and Exhibition, SPE30595MS. https://doi.org/10.2118/30595MS Brown RJ, Korringa J (1975) On the dependence of elastic properties of a porous rock on the compressibility of the pore fluid. Geophysics 40 (10): 608–616. https://doi.org/10.1190/1.1440551 Digby PJ (1981) The effective elastic moduli of porous granular rocks. Journal of Applied Mechanics 48: 803–808. https://doi.org/10.1115/1.3157738 Dvorkin JP, Nur AM (1996) Elasticity of highporous sandstones: Theory for two North Sea data sets. Geophysics 61(5): 1363–1370. https://doi.org/10.1190/1.1444059 Eshelby JD (1957) The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proceedings of Royal Society London 241(1226): 376–396. https://doi.org/10.1098/rspa.1957.0133 Graebner M (1992) Planewave reflection and transmission coefficients for a transversely isotropic solid (short note). Geophysics 57(11): 1512–1519. https://doi.org/10.1190/1.1443219 Hashin Z, Shtrikman S (1962) A variational approach to the theory of effective magnetic permeability of multiphase materials. Applied Physics 33: 3125–3131. https://doi.org/10.1063/1.1728579 Hashin Z, Shtrikman S (1963) A variational approach to the elastic behavior of multiphase materials. Mechanics and Physics of Solids 11: 127–140. https://doi.org/10.1016/00225096(63)900607 Hill R (1963) Elastic properties of reinforces solids: some theoretical principles. Mechanics and Physics of Solids 11: 357–372. https://doi.org/10.1016/00225096(63)90036X Hornby BE, Schwartz M, Hudson JA (1994) Anisotropic effectivemedium modeling of the elastic properties of shales. Geophysics 59(10): 1570–1583. https://doi.org/10.1190/1.1443546 Hsu K, Esmersoy C, Schoenberg M (1988) Seismic velocities and anisotropy from highresolution sonic logs. SEG Technical Program Expanded Abstracts 7(1): 1359. https://doi.org/10.1190/1.1892475 Hudson JA (1980) Overall properties of a cracked solid. Mathematical Proceedings of the Cambridge Philosophical Society 88: 371–384. https://doi.org/10.1017/S0305004100057674 Hudson JA (1981) Wave speeds and attenuation of elastic waves in material containing cracks. Geophysical Journal International 64: 133–150. https://doi.org/10.1111/j.1365246X.1981.tb02662.x Imhof MG (2003) Scale dependence of reflection and transmission coefficients. Geophysics 68(1): 322–336. https://doi.org/10.1190/1.1543218 Ismail A, Ewida HF, AlIbiary MG, Zollo A (2020) Application of AVO attributes for gas channels identification. Petroleum Research 5(2): 112–123. https://doi.org/10.1016/j.ptlrs.2020.01.003 Jenkins J, Johnson D, La Ragione L, Makse H (2005) Fluctuations and the effective moduli of an isotropic, random aggregate of identical, frictionless spheres. Journal of the Mechanics and Physics of Solids 53: 197–225. https://doi.org/10.1016/j.jmps.2004.06.002 Kinoshita N, Mura T (1971) Elastic fields of inclusions in anisotropic media. Physica Status Solidi 5(3): 759–768. https://doi.org/10.1002/pssa.2210050332 Liner CL (2006) Layerinduced seismic anisotropy from full wave sonic logs. SEG Technical Program Expanded Abstracts 25: 159–163. https://doi.org/10.1190/1.2369826 Li L, Ma J, Wang H, Tan M, Cui S, Zhang Y, Qu Z (2017) Study of shear wave velocity prediction during CO_{2}EOR and sequestration in Gao 89 area of Shengli Oilfield. Applied Geophysics 14(3): 372–380. https://doi.org/10.1016/10.1007/s1177001706385 Ma J, Morozov IB (2010) AVO modeling of pressuresaturation effects in Weyburn CO_{2} sequestration. The Leading Edge 29 (2): 178–183. https://doi.org/10.1190/1.3304821 Mavko G, Mukerji T, Dvorkin J (2020) The rock physics handbook. Cambridge University Press, Cambridge, 222–233 Mindlin RD (1949) Compliance of elastic bodies in contact. J. Appl. Mech 16(3): 259–268. https://doi.org/10.1115/1.4009973 Mura T (1991) Micromechanics of defects in solids. Kluwer Academic Publisher, Amsterdam, The Netherlands, 368–386 Norris AN, Johnson DL (1997) Nonlinear elasticity of granular media. ASME Journal of Applied Mechanics 64: 39–49. https://doi.org/10.1115/1.2787292 Sassen R, Losh LS, Cathles L, Roberts H, Whelan J, Milkov A, Sweet S, DeFreitas D (2001) Massive veinfilling gas hydrate: Relation to ongoing gas migration from the deep subsurface in the Gulf of Mexico, Marine and Petroleum. Geology 18(5): 551–560. https://doi.org/10.1016/S02648172(01)000149 Schoenberg M, Protázio J (1992) "Zoeppritz" rationalized and generalized to anisotropy. Journal of Seismic Exploration 1: 125–144. https://doi.org/10.1121/l.2029011 Smith R, Bakulin A, Jervis M, Hemyari E, Alramadhan A, Erickson K (2018) 4D seismic monitoring of a CO_{2}EOR demonstration project in a desert environment: acquisition, processing and initial results. SPE Kingdom of Saudi Arabia Annual Technical Symposium, 192311MS. https://doi.org/10.2118/192311MS Sondergeld CH, Rai CS (2011) Elastic anisotropy of shales. Leading Edge 30 (3): 324–331. https://doi.org/10.1190/1.3567264 Papageorgiou G, Amalokwu K, Chapman M (2016) Theoretical derivation of a Brielike fluid mixing law. Geophysical Prospecting 64(4): 1048–1053. https://doi.org/10.1111/13652478.12380 Patil Pramod D, AlQasim Abdulaziz S, AlZayani Alia I, Kokal Sunil L (2023) Endtoend surface and subsurface monitoring and surveillance for the onshore and offshore CCS or CCUS Projects. Paper Presented at the Offshore Technology Conference, Houston, Texas, USA. https://doi.org/10.4043/32172MS Qian K (2017) Anisotropic rock physics modeling and brittleness analysis of shale reservoir. Doctoral thesis, China University of Petroleum, Beijing, 15–18 Rüger A (1996) Reflection coefficients and azimuthal AVO analysis in anisotropic media. Colorado School of Mines, Denver, The United States, 65–82 Rüger A (1997) Pwave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry. Geophysics 62(3): 713–722. https://doi.org/10.1190/1.1444181 Thomsen L (1986) Weak elastic anisotropy. Geophysics 51(10): 1954–1966. https://doi.org/10.1190/1.1442051 Wang H, Ma J, Li L, Jia L, Tan M, Cui S, Zhang Y, Qu Z (2017) Timelapse seismic analysis for Gao89 area of CO_{2}EOR project in SINOPEC Shengli Oilfield, China. Energy Procedia 114: 3980–3988. https://doi.org/10.1016/j.egypro.2017.03.1530 Wang P, Chen X, Li J, Wang B (2019) Accurate porosity prediction for tight sandstone reservoir: a case study from north China. Geophysics 85(2): 1–71. https://doi.org/10.1190/geo20180852.1 Wang P, Li J, Chen X, Wang B (2020) Joint probabilistic fluid discrimination of tight sandstone reservoirs based on Bayes discriminant and deterministic rock physics modeling. Journal of Petroleum Science and Engineering 191: 107218. https://doi.org/10.1016/j.petrol.2020.107218 Wang P, Cui YA, Liu J (2022) Fluid discrimination based on inclusionbased method for tight sandstone reservoirs. Surveys in Geophysics 43(5): 1469–1496. https://doi.org/10.1007/s10712022097125 Walton K (1987) The effective elastic moduli of a random packing of spheres. Journal of the Mechanics and Physics of Solids 35: 213–226. https://doi.org/10.1016/00225096(87)900366 White DJ (2009) Monitoring CO_{2} storage during EOR at the WeyburnMidale Field. The Leading Edge 28(7): 838–842. https://doi.org/10.1190/1.3167786 White DJ (2013) Seismic characterization and timelapse imaging during seven years of CO_{2} flood in the Weyburn field, Saskatchewan, Canada International Journal of Greenhouse Gas Control 16(11): S78–S94. https://doi.org/10.1016/j.ijggc.2013.02.006 Wood AW (1955) A text of sound. The MacMillan Co, New York, 236–258 Wu F, Li J, Geng W, Tang W (2022) A VTI anisotropic media inversion method based on the exact reflection coefficient equation. Frontiers in Physics 10: 926636. https://doi.org/10.3389/fphy.2022.926636 Yuan H (2007) Study of anisotropic rock physics model and application. China University of Geosciences, Beijing, 56–49 Zhao X (2017) Methodologies of prestack seismic inversion for shale gas play. China University of Petroleum, Beijing, 29–43 Zoeppritz K (1919) On the reflection and propagation of seismic waves. Gottinger Nachrichten 1: 66–84