Chapter 12 Seismic Reflectivity in Carbon Dioxide Accumulations : A Review

It is widely recognized that the continuous increase in the concentration of carbon dioxide (CO2) in our atmosphere is a major cause of global climate change [1]. Consequently, in accordance with the objectives of the Kyoto Protocol, many carbon dioxide capture and storage projects are being developed worldwide to reduce and stabilize the emission of this greenhouse gas into the atmosphere. The geological storage of CO2 in many cases becomes a feasible option to accomplish this goal, giving rise to the science of carbon sequestration, a challenging task for governments, scientists and engineers [2].


Introduction
It is widely recognized that the continuous increase in the concentration of carbon dioxide (CO 2 ) in our atmosphere is a major cause of global climate change [1]. Consequently, in accordance with the objectives of the Kyoto Protocol, many carbon dioxide capture and storage projects are being developed worldwide to reduce and stabilize the emission of this greenhouse gas into the atmosphere. The geological storage of CO 2 in many cases becomes a feasible option to accomplish this goal, giving rise to the science of carbon sequestration, a challenging task for governments, scientists and engineers [2].
The main targets for geological disposal of CO 2 are depleted hydrocarbon reservoirs and deep saline aquifers. While the latter are more numerous, their characterization requires detailed studies because they are typically not explored for prospecting. The geological sequestration of carbon dioxide requires careful prior study and subsequent monitoring to prevent this gas from leaking to the surface. In general, the significant contrast between the physical properties of natural reservoir fluids and those of carbon dioxide allows the utilization of time lapse seismic methods to monitor the evolution of the injected CO 2 volume.
While it is accepted that time lapse surveys are able to monitor the presence or absence of CO 2 within a geological formation, their ability to quantify the saturation, state and volume of this fluid within the reservoir still needs research efforts from the geophysical community. This makes it necessary to search for reliable seismic indicators. Theoretical and numerical modeling provides, in this sense, a tool to explore useful correlations between the seismic response and relevant parameters of the CO 2 repository.
The significance of the seismic reflectivity has long been recognized by many authors, resulting in the publication of numerous works for different constitutive models and applications. In close connection with this, the analysis of seismic amplitude variation with angle, hereafter called AVA is frequently used in reservoir geophysics to obtain information about the rocks and pore fluids. The behavior of the reservoir reflectivity for different overall CO 2 saturations, distribution types, layer thicknesses, frequencies and CO 2 physical states controls the amplitude of seismic wave reflections and strongly conditions the detectability of the injected CO 2 volume, as reported by [3], [4] and [5] among others.
With this idea, using rock physics models, in this chapter we focus on the theoretical analysis of the properties and variations of the seismic reflectivity and related parameters in a porous rock partially saturated with CO 2 and brine. In our tests we model the seismic effects of different CO 2 accumulations formed below impermeable seals.
We calibrate our models using information on the Utsira formation, a shallow saline aquifer in the Sleipner field (offshore Norway). Millions of tonnes of CO 2 have been separated since October 1996 from natural gas and re-injected into the aquifer, consisting of a highly porous and unconsolidated sandstone, with several thin intra-reservoir shale layers. These shale intervals act as temporary seals causing accumulations of CO 2 beneath them [6], [7], whose saturations increase with time. A 50-100 m thick high velocity layer overlies the sand, forming the reservoir caprock, known as the Nørland formation. It is mainly composed of shales, siltstones and mudstones [8]. Another high-velocity shale layer of about 5 m thick is located about 10 m below the top of the sand.
As pointed out by [9] the topmost CO 2 layer is of special interest for monitoring the Sleipner injection site, since its growth reveals the total upward flux and its changes over time. Thus, the analysis of the reflectivity associated with this kind of CO 2 accumulation is an important topic in CO 2 sequestration problems. The present book chapter focuses on this subject and is structured in three main sections considering different aspects of the problem.
First, in Section 1 we present a description of the models and rock physics tools used in the subsequent sections. Next, in Section 2 we model the compressional P-wave reflection coefficient at the interface between a caprock and a permeable porous layer partially saturated by a mixture of brine and CO 2 under liquid, supercritical and gaseous conditions. For this analysis we consider a simple model consisting of two halfspaces, combining a fluid substitution procedure and a Gassmann-Hill formulation to take into account variable spatial distributions of the fluids. We perform a sensitivity analysis of the standard AVA coefficients (intercept, gradient and curvature) in the near offset range, to investigate their correlations with CO 2 saturation and its thermodynamic state within the geologic formation. Also, we analyze the effect of modeling the CO 2 by means of simple and complex equations of state.
The two-halfspace representation considered in Section 2 is valid when the accumulation thicknesses are larger than the involved seismic wavelengths. When this assumption is not valid, the interference between the multiple waves reflected at the different boundaries give rise to strong frequency dependence, reflectivity dispersion and tuning effects which must be taken into account. To do this, in Section 3 we compute the generalized P-wave reflection coefficients by means of a reflectivity method. This leads us to the study of frequency dependent AVA and amplitude vs. frequency AVF, which are topics of great interest due to the increasing use of spectral decomposition techniques on seismic data. This also allows us to analyze the pattern of maxima in the reflectivity amplitude and their associated peak frequencies, which are strongly related to the thickness of the CO 2 layer.
Seismic monitoring of geological CO 2 injection sites is mostly based on the seismic reflections coming from high saturation CO 2 accumulations. This is due to their large seismic amplitudes and good signal to noise ratio. However, low-saturation zones with dispersed CO 2 , or saturation transition zones may have an important role in the propagation of waves within the reservoir, giving rise to amplitude and phase changes of the seismic signals. Transition zones have been studied by different authors in the geophysical literature considering linear velocity trends with depth and constant density [10], [11]. Therefore, using the parameters of the Sleipner field, in Section 3.3 we model the reflectivity response of a CO 2 transition layer, defined by a given linear vertical CO 2 saturation profile, which results in a non-linear velocity trend with depth.
In general, emphasis is placed on establishing correlations between the seismic reflectivity and related attributes (intercept, gradient, curvature, peak frequencies) and the overall CO 2 saturation, its physical state and thickness of the accumulation. These results are intended to help in the understanding of the expectable variations in a seismic time lapse study. They can be extended to other CO 2 repositories with proper calibration of the rock and fluid properties.
Many of the results and procedures presented in this chapter are a revision and extension of those published by us in [4], [12], [13], [14] and [15].

Theoretical framework and assumptions
In this section we summarize the modeling tools used later for the applications. First we explain the procedure for the calculation of the elastic properties of the partially saturated rocks for variable CO 2 saturation. Next we describe the laws for the computation of the physical properties of the CO 2 and brine for the different temperature and pressure states. Finally, we outline the reflectivity method for the computation of the generalized seismic reflectivity of the layered model

Elastic properties of CO 2 bearing rocks
To begin with the description of the model we consider that after the injection, a volume of carbon dioxide occupies part of the pore volume of a geologic reservoir which at the pre-injection state was fully saturated with brine. For simplicity we assume that CO 2 displaces, without dissolution, part of the in-situ brine giving rise to a two-phase fluid saturation. From now on the corresponding saturations of free CO 2 and brine are denoted as Sg and S br , respectively, so that S br + Sg = 1.
For the study of partially saturated rocks it is necessary to compute the bulk density and elastic coefficients (bulk and shear modulus) of the fluid saturated medium. The bulk density ρ is given by where φ is the rock porosity, ρs is the mineral grain density and ρg, ρ br are CO 2 and brine densities, respectively.
The formulation and solution of the energy and amplitude splitting problem when a monochromatic plane compressional wave strikes obliquely at a plane interface between two porous saturated media has been discussed by various authors, such as [16], [17] and [18]. In those papers the mechanical behavior of the porous rock was described using the classic constitutive relations and equations of motion of Biot [19], [20]. In the following sections our aim is to study the reflection of elastic waves and related parameters, for frequencies f within the common seismic range f 120Hz. When a low frequency seismic wave propagates through a fluid saturated rock, due to fluid viscosity the solid and the fluid move in phase, so the medium behaves as an effective medium.
One of the most significant models to estimate the effective bulk modulus and the seismic velocity of a porous fluid saturated rock is to use Gassmann's relations [21] which can be written in the form: In this equation the mechanical behavior of the fluid saturated porous medium is assumed to be elastic and isotropic. The coefficients Km, Ks and K f are the bulk modulus of the dry matrix, the mineral grains and the pore fluid, respectively. The physical properties of the two-phase brine-CO 2 fluid are computed using an effective fluid whose density is given by their weighted average and its compressibility by the isostress Reuss average of individual fluid bulk moduli [22]: where K br and Kg are the bulk moduli of brine and CO 2 . In these equations, it is assumed that the mixture of CO 2 and brine at the pore scale can be treated as a viscous single phase effective fluid. However, as pointed out in [23], this approach is strictly valid only when the pore fluids are uniformly mixed at very small scales so that the different wave-induced pressure increments in each fluid have time to diffuse and equilibrate during a seismic period. This critical distance is the so-called diffusion length. The presence of fluid heterogeneities over scales greater than this length (in which wave induced pore pressure gradients cannot equilibrate quickly), gives rise to a patchy saturation. It has been shown that for low frequencies such as those used in seismic exploration the effective modulus of a rock with patches of brine and CO 2 of arbitrary geometry is given by [24], [23] where K G br and K G g are the Gassmann's moduli (2) of the rock fully saturated with brine and CO 2 . The coefficient µ denotes the shear modulus of the rock, equal to that of the rock matrix, since the rigidity of a rock does not change due to the saturant fluid.
The existence of heterogeneous fluid distribution may give rise to mesoscopic wave attenuation and velocity dispersion phenomena in the seismic frequency range. There are many studies about numerical modeling of these effects, e.g. [25], [26], [27], showing that they are strongly dependent on the shapes and characteristic lengths of the patches. Given that these parameters are rarely known, to avoid dealing with this uncertainty, in our analysis we assume elastic layers where no attenuation-dispersion phenomena, associated with the patchy fluid distribution take place. However, these effects can be included in our formulation using appropriate viscoelastic models in the frequency domain. Because of the classic correspondence principle, this implies the replacement of the real constant elastic moduli by complex frequency dependent moduli [20].
The patchy (Hill) and uniform (Gassmann) fluid distributions are respectively upper and lower bounds for the seismic wave velocities [23]. Taking into account the uncertainties in the knowledge of the in-situ CO 2 -brine distribution, in our models we assume that for each saturation it is reasonable to compute the average between both velocities.
It is important to mention the assumption of no chemical interactions between the pore fluids and the frame, which allows us to employ the fluid substitution procedure to consider that the pore space is saturated by brine and CO 2 in variable proportions. Also, since the amount of CO 2 dissolved in brine is a negligible fraction [28], [29], we do not take this effect into account. Nevertheless, when CO 2 is injected in oil reservoirs this effect can not be neglected and should be included in the physical models, as shown in [4] and [28].
As pointed out in [8], no systematic variations in fluid pressure are observed in the Sleipner field. Thus, variations of rock elastic properties with effective pressure (related to the difference between confining and pore pressure), are not taken into account in the model. However, using appropriate effective pressure laws these effects can be included in the computations, as explained in [18].

Physical properties of carbon dioxide and brine
Depending on the in-situ pressure and temperature conditions the CO 2 can exist in the subsurface in different phases. We recall that the critical point for CO 2 occurs at a temperature Tc = 31.1 • C and a pressure Pc = 7.39 MPa. For temperatures higher than Tc and pressures higher than Pc the carbon dioxide is in a supercritical state, where it is compressible like a gas but with the density of a liquid. This characteristic of CO 2 is particularly relevant for its underground storage since supercritical CO 2 can fill the available pore volume with minimum buoyancy effects [30]. Temperatures and pressures near the critical point commonly occur in applications involving CO 2 , such as enhanced oil recovery techniques and sequestration projects [31]. However, as pointed out by [30], the depth at which CO 2 supercritical conditions are present is highly variable and strongly dependent on surface temperature and geothermal gradients, even within a single basin. In addition, the pressure regime of the basin (normal or abnormal), is also very important and is related to its geologic history, existence of sealing faults, permeability barriers and the occurrence of overpressure generation mechanisms [32].
For the examples, the density and bulk moduli of brine for given in-situ temperature and pressure conditions are computed using the semi-empirical relations proposed by Batzle and Wang (BW) [29]. For the computations we consider a typical brine salinity of 50000 ppm. The corresponding properties of CO 2 can be computed using some of the many equations of state (EoS) developed for real gases such as the two-parameter van der Waals (vW) equation [35] and the Peng and Robinson (PR) equation [34], or some more specific one such as Duan, Moller and Weare (DMW) [33], which involves fifteen parameters. To our knowledge, there is no full agreement in the geophysical literature about which EoS is the most appropriate to represent CO 2 properties under the thermodynamic conditions found in geologic reservoirs. Therefore, we use different models to analyze the effect of the EoS in the seismic magnitudes of interest. With that purpose we selected three different temperature Table 1. Physical properties of CO2 and brine under different thermodynamic conditions. The density and bulk moduli of CO2 are obtained using three equations of state: DMW [33], PR [34] and vW [35]. The CO2 bulk moduli estimations are corrected for adiabatic conditions. The properties of brine are computed using BW laws [29] for a salinity of 50000 ppm.  Table 1 and will be used in the following sections. The temperature and pressure of the supercritical state are in agreement with the conditions reported in [8] for Utsira sandstone.
It must be pointed out that the different CO 2 bulk moduli estimations in Table 1, to be used for the computations, were corrected to represent adiabatic conditions. The reason is that the fluid compression due to the passage of a wave is fast, so that the process is not isothermal [29], [37].

Calibration of the rock physics model
To calibrate our elastic model for the Utsira sandstone, we used the measured bulk density and the compressional and shear wave velocities in the preinjection state (i.e. under full brine saturation), given in [8] and [38]: ρ b = 2050 kg/m 3 , Vp =2050 m/s and Vs = 640 m/s. Using the parameters given for the pore water ρ br = 1040 kg/m 3 , K br = 2.305 GPa and the average sandstone porosity φ = 0.37, the mineral grain density results in ρs = 2643 kg/m 3 . We also performed a Gassmann-inverse calculation to obtain: Km = 2.569 GPa and µ = 0.84 GPa. For the shale layers we obtained the elastic parameters using the reported values: Vp = 2270 m/s, Vs = 850 m/s and ρ b = 2100 kg/m 3 .

The generalized reflectivity of a layered medium
The computation of the reflection and transmission coefficients of elastic waves propagating through layered media is a subject of interest in different fields such as seismology, seismic prospecting and underwater acoustics. The first numerical procedures date back to the pioneering works of Thomson [39] and Haskell [40], who proposed matrix methods that transfer stresses and displacements through successive layers. Since then, different approaches and applications have been presented.
In this work we perform an implementation of the reflectivity method [41], [42], which is based on the continuity of particle displacements and stress components at the interfaces between sets of plane layers embedded between two halfspaces. This method builds up the reflection and transmission matrices iteratively by starting at the top of a lower bounding halfspace and adding one layer per iteration until the total stack response is constructed. This procedure is intuitively simple and exact [43], taking into account all internal reverberations. The recursion algorithm is were R j and T j are the total-reflection and total-transmission matrices of interface j.
are the upward and downward layer reflection and transmission matrices respectively. These matrices depend on the bulk density, P-and SV-wave seismic velocities and ray angles from interface j. E j is the layer phase-shift diagonal matrix and (I − R u j−1Rj ) −1 is known as the reverberation operator [44], with I being the identity matrix.
The recursion starts at the base of the layering j = n with Rn = R d n = 0 and Tn = T d n = I and continues to j = 1. When the top of the layering is reached, the generalized P-wave reflection coefficient for scaled displacements is obtained from Rpp(f , θ) = R 0 [1,1] . This coefficient is a complex function of frequency f and incidence angle θ.
For a complete description of the matrices involved in this scheme the reader is referred to [44].

AVA analysis at the top of a CO 2 accumulation
In this section we model the seismic compressional wave reflection coefficient and related AVA parameters at the interface between Utsira sandstone, partially saturated by mixtures of CO 2 and brine, overlain by a shale caprock, with the properties described in Section 1.3.
Following the classic approach, for the present applications both media are considered as elastic halfspaces, an assumption valid for layers of thicknesses larger than the wavelengths of the incident waves. Under these conditions the reflection coefficient are real and frequency independent. In this particular configuration, our reflectivity method reproduces the classic results of Zoeppritz [45] for two halfspaces.
In the near offset domain θ < 60 • or below the critical angle, we can assume that the Rpp reflection coefficient as a function of the incidence angle θ can be approximated using the expression of Shuey [46]: where A, B and C are known as the AVA coefficients. A is called the intercept, B the gradient and C the curvature [47], [48]. The intercept is equal to the normal incidence reflection coefficient and is controlled by the contrast in acoustic impedance between both media. The gradient is related to the contrast in density and in compressional P and shear SV wave  velocities [49]. The curvature parameter is important at far offsets and near critical angles [44]. Equation (6) is used in this section to carry out a parametric analysis of the AVA coefficients or to study their sensitivity by implementing a standard fitting procedure on the results obtained for Rpp(θ).
In Figure 1 we plot the compressional elastic wave velocity for the partially saturated CO 2 -brine Utsira sandstone in the low frequency range. The parameters of the CO 2 at the liquid, supercritical and gaseous states correspond to the PR EoS and are listed in Table 1, along with those of the brine. As mentioned in Section 1, our velocity model for the sandstone corresponds to the Gassmann-Hill average for each saturation state. It is important to note the decreasing tendency of the wave velocity for increasing CO 2 saturation, which becomes less marked as saturation increases, a feature that will also be observed in the AVA coefficients. The velocity model for the supercritical CO 2 is in close agreement with that presented in [8] and [38] for the Utsira sandstone. According to this model, the seismic velocity can discriminate the liquid from the other two physical states, particularly for saturations lower than about 40%.
Next, in Figure 2 a)-d) we show the compressional reflection coefficient versus incidence angle from 0 • to 50 • , for fixed values of CO 2 saturation in the range 0%-100%. The temperature is 36 • C and pressure 10 MPa, corresponding to the supercritical conditions given in Table 1, for the PR equation. As expected, an abrupt change in the reflectivity from the pre-injection to the post-injection state is observed, due to the high contrast between the physical properties of CO 2 and brine in this state. As expected, the CO 2 saturation decreases the impedance of the lower medium with respect to that of the upper one, resulting in a negative reflection coefficient for normal incidence, becoming more negative for higher angles and higher saturations. These curves are consistent with those published in [50]. The curves corresponding to liquid and gaseous states show very similar behavior, and consequently they are not shown. This means that the discrimination of the physical state from the AVA behavior may be difficult. This subject will be analyzed in more detail in terms of the AVA coefficients. According to the well known AVA classification [48], [51], the curves for the preinjection case and for Sg = 10% (in Figure 2 b) correspond to the AVA class IV. The set of AVA curves for higher saturations shows a very flat tendency for near angles (i.e. with very low gradients B), which can be observed in more detail in Figures 2 c) and d). Its classification will be discussed later in Figure 4.
To analyze the sensitivity of the AVA parameters we plot in Figure 3 the predicted variations of A, B and C with increasing saturation for the liquid, supercritical and gaseous conditions. We observe that the intercept A is negative, decreasing and monotone for the three physical states under consideration, being a sensitive attribute throughout the saturation range. The gradient B shows a trend similar to that of the intercept except for the gaseous state, in which stabilization is observed. This can be useful to distinguish this particular state from the other two.
The curvature attribute C, is the least sensitive of the three, both to physical state and saturation level. We only note appreciable changes for liquid CO 2 and for saturations in the range from 0 to about 30%. However, the correct determination of such changes may be strongly limited by the seismic resolution. In connection with this, Brown et al. [3] pointed out that AVA variations on the order of 5% are seismically detectable. Taking that numerical threshold into account we can state that, for this model, variations in the parameter C do not contain much information with respect to variations in the CO 2 content. Consequently, an AVO time-lapse analysis can be based solely on the coefficients A and B, where A is the most sensitive for the determination of the CO 2 saturation state. To inspect the behavior of the reflectivity in terms of the AVA classification, in Figure 4 we present typical intercept-gradient crossplots (B vs. A) to observe their trends for variable saturations and physical states. According to our model, while for gaseous CO 2 the AVA class is IV in all the saturation range, under liquid and supercritical conditions a change to class III can be expected for very high saturations. We also observe that since from intermediate to high gas saturations, the B − A trajectories of each CO 2 state are different, a time-lapse crossplot analysis may bring useful information about the physical state of the underground CO 2 . In order to analyze the influence of the state equations in our model estimates, in Figure 5 we display the coefficients A and B as a function of CO 2 saturation, for van der Waals [35], Peng & Robinson [34] and Duan et al. [33] EoS for the three physical states in Table 1. For gaseous and supercritical states the coefficients are very similar in all the saturation range. The most important discrepancies are observed when CO 2 is at liquid state. From these experiments we conclude that for seismic modeling applications the choice of the EoS will not significantly affect the results.
The sensitivity and robustness of the AVA coefficients to the presence of low fractions of methane (CH 4 ), which is a frequent impurity in CO 2 injection sites, was considered in [15].

Finite thickness CO 2 accumulations: frequency dependent AVA
So far we have modeled the CO 2 reflectivity response in terms of the simple single layer representation. However, due to the spatial variations in petrophysical properties and natural permeability barriers, the CO 2 tends to form accumulations of finite thickness.
In exploration geophysics, it is generally accepted that when the thickness of a layer is smaller than the predominant seismic wavelengths, the layer is said to be thin. Also, a common measure for the vertical seismic resolution is given by a quarter of the dominant wavelength [53]. As explained in [6], [8] and [54], the CO 2 accumulations within the Utsira sandstone in the Sleipner field can be regarded as thin layers. The interference effects between the multiple waves reflected at the different boundaries give rise to strong frequency dependence in the reflectivity of the medium, a phenomenon known as reflectivity dispersion [11], [55]. In this case the dispersive character has a geometric origin and is due to layer reverberation; however reflectivity dispersion can also take place at interfaces between highly dissipative media. The constructive interference between some frequency components of the up-going and down-going reflections at the physical boundaries gives rise to tuning effects. In the context of CO 2 sequestration these effects have been analyzed in several papers, such as [5], [8], [9], and [38].
As pointed out in [11] and [55], the modern use of time-frequency analysis and spectral decomposition tools [56] has shown that reflection events in practice are always frequency dependent. The spectral decomposition techniques allow obtaining and analyzing the seismic response of a medium at a given frequency. This makes it possible to investigate the reflected amplitude variation with angle at different frequencies (spectral AVA) and the amplitude variation with frequency (AVF ) at a fixed angle.
At the time of this writing and to our knowledge, the AVF characteristics of CO 2 layers have not been modeled yet. With this motivation, in this section we investigate the existing correlations between reflectivity, frequency, saturation, physical state and thickness for a simple accumulation model. From these experiments our aim is to assess whether time lapse changes in the reflected amplitudes can be correlated with changes in such parameters, which are of utmost importance in CO 2 monitoring.

Description of the model
Here we present another simple but useful approach to model the reflectivity of the topmost CO 2 layer in the Sleipner field. For the computations we use the sandstone and caprock parameters described in Section 1.3. To simulate the supercritical conditions of the CO 2 injected into Utsira sandstone, carbon dioxide density and bulk moduli are calculated using the PR EoS at a temperature of 36 • C and pressure of 10 MPa (Table 1).
According to the scheme used for the reflectivity method (described in Section 1.4), the upper halfspace plays the role of the caprock and the lower halfspace is the fully brine saturated Utsira sandstone. Below the caprock, the injected carbon dioxide is assumed to form a layer in the sandstone of thickness h with spatially constant CO 2 saturation Sg. The case of vertical variations in the CO 2 saturation or saturation transitions will be considered in Section 3.3. Thus for this model the generalized reflection coefficient is a complex, frequency-dependent function of the incidence angle, the CO 2 saturation and thickness h in the CO 2 bearing layer.

Behavior of reflectivity with saturation, frequency and thickness
In Figure 6 we display the generalized reflection coefficient Rpp vs. incidence angle for a frequency of f = 50 Hz, a saturation of Sg = 100 % and thicknesses of h = 0, 5, 10 m, considering also the limit h → ∞. For this CO 2 saturation the P-wave velocity (see Figure  1) is Vp ≈ 1410 m/s, so that the associated wavelength λ = Vp/f is about 28 m and the vertical resolution limit, λ/4, at this frequency results in about 7 m. So, in these examples the 5 m layer is ultra-thin and the 10 m layer is thin.
When h = 0 and h = ∞, the model reduces to the classic two halfspace model, with a real reflection coefficient (shown in Figures 6 a) and b)), which can be analyzed in terms of the AVA classification. As can be seen in Figures 6 b) and c), the spectral AVA response for finite thickness layers differs considerably from that of the halfspaces, as expected. For 0 < h < ∞ the AVA character is highly dependent on frequency, a dependence that can be clearly observed in Figure 7. In particular, the curves for the thicknesses h = 5 m and h = 10 m show a periodic behavior; meaning that the different frequency components of a seismic signal are reflected with different amplitudes. To understand this periodicity on a theoretical basis, it is useful to consider the analytical result for the reflectivity of an acoustic layer embedded between two halfspaces, treated in [52] and [53]. From that solution it can be shown that the function Rpp(f , θ) has multiple maxima or peaks at frequencies f n max and minima or notches located at frequencies f n min , given respectively by [15] f n max = The occurrences of these extrema are directly proportional to the compressional velocity within the layer Vp (shown in Figure 1) and inversely proportional to its thickness h, where θ L is the ray angle inside the layer. The frequencies at which the reflected amplitude has maxima are known as tuning frequencies [58] or peak frequencies.
The first peak frequency f 0 max for normal incidence, hereafter denoted as fp, results in Equation (8) is a very useful relation for thickness determination by means of the fp attribute, as shown in [38] in the context of CO 2 sequestration.
For an arbitrary number of layers with different physical properties and variable thicknesses the generalized reflection coefficient also displays a pattern of maxima and minima. However, for such cases there are no simple analytical expressions and consequently a numerical procedure such as the one presented in this work enables predicting the spectral characteristics of their reflectivity.
For the particular cases h = 0 and h → ∞, plotted in Figures 7 a) and d), the reflectivity for normal incidence is real, constant and related to the impedance contrast between the caprock and the fully brine saturated sandstone and the fully CO 2 saturated layer. In Figures 7 b) and c), the location of the peak close to 70 Hz for h = 5 m and the peaks close to 35 Hz and 105 Hz for h = 10 m can be verified. Note also that for the continuous component f = 0 Hz, these AVF curves converge to the response in Figure 7 a), given that for such a long wavelength the layers h = 5 m and h = 10 m are invisible to the seismic waves.  The modulus of Rpp versus frequency and thickness for normal incidence is shown in Figure  8 a), in which the reciprocal relation between peak frequencies and thicknesses is clearly displayed. As seen, at each thickness we have a different peak frequency. To study this dependency in more detail, in Figure 8 b) we display the first peak frequency versus thickness, in which it is instructive to note that for a frequency of 35 Hz the corresponding thickness h is about 10 m. This is the case studied in [38], where spectral decomposition is used at Sleipner to map the topmost CO 2 accumulation. For very small thicknesses Figure 8 b) also indicates that the peak frequencies are well above the standard seismic frequency range, so high-resolution seismic acquisition is essential for these layers.
To analyze the sensitivity of the previous results to CO 2 saturation, next in Figure 9, we plot the modulus of Rpp versus frequency and saturation for thickness h = 10 m and incidence angles 0 • , 30 • and 60 • . The peaks of the reflectivity are located in the high CO 2 saturation range, since in that case the impedance contrast between the caprock and the sandstone is maximum. The sensitivity to Sg depends on frequency showing bands of very weak reflectivity (in blue scale) and higher amplitudes (in red scale). The high reflectivity zones, which at the same time are sensitive to saturation, are located around the first and second peak frequencies for this example (h = 10 m), i.e. near 35 Hz and 105 Hz. These conclusions are valid for normal incidence θ = 0 • and also for larger angles, although a slight shift of the peak frequencies to higher values is observed, as expected from equation (7).

The influence of saturation transitions
The seismic monitoring of carbon dioxide in geologic reservoirs is mostly focused on the characterization of accumulations of high saturation, due to their large seismic amplitudes. Nevertheless, low-saturation zones with dispersed CO 2 , or saturation transitions may have an important role in the propagation of waves within the reservoir, giving rise to amplitude and phase changes of the seismic signals. Our aim is to analyze whether the presence of a saturation transition below the main accumulation substantially alters the previous modeling results. At the same time we investigate the influence of the vertical extent of the transition zone on the reflectivity of the reservoir.
To study these effects, we consider a modification of the model used in Section 3.2, by introducing a variable CO 2 saturation-depth profile. This model supposes that the migration of the injected CO 2 is controlled by buoyancy and that the reservoir permeability is homogeneous and isotropic, resulting in a laterally constant gas saturation field. The saturation-depth profile is inspired by the axisymmetric flow simulations presented in [9].
For the following examples we represent the saturation Sg(z) with a piecewise linear function of the vertical depth z. The total thickness of the CO 2 bearing layer is given by where h 0 is the thickness of a constant saturation zone below the caprock, with Sg = 100 % and h t being the thickness of the saturation transition. The vertical variation in saturation Sg(z) in turn produces a linear density profile ρ(z) and a non-linear velocity with depth Vp(z), as shown in Figure 10. The physical conditions of the CO 2 are those described in Section 3.1. To compute the generalized reflection coefficient using the reflectivity method, Figure 11. Normal incidence AVF for transition thicknesses ht = 0, 5, 10 and 20 m.
the upper halfspace of the model is the caprock, followed by the constant saturation layer with h 0 = 10 m, the saturation transition with variable h t and a fully brine saturated lower halfspace. The transition is discretized in a finite number of layers of 0.5 m with decreasing CO 2 saturation with depth.
In Figure 11 a)-d) we compare the normal incidence AVF response for transitions of different thicknesses given by h t /h 0 = 0, 0.5, 1 and 2, that is h t = 0, 5, 10 and 20 m. The overall bulk saturation within the transition, given bȳ is kept constant atS = 50 % in the four cases. Due to the layered geometry and variable saturation of the model, the reflectivity vs. frequency shows a pattern of maxima and minima, with decreasing peak amplitudes with frequency. The presence of the saturation transition shifts the first peak to lower frequencies within the seismic band. In general for higher thicknesses h t we observe lower values of the first peak frequency, which indicates that the inverse proportionality between these parameters still holds. Although the theoretical relation in equation (7) is not applicable to this multilayer geometry, the distinct periodicities can help to discriminate thin and thick transition zones.
The effect of the saturation transition can be noted in more detail in Figure 12, where we plot the first peak frequency fp vs. thickness h for the case with and without saturation transition below the main accumulation. The discrepancies between the curves indicate that the CO 2 saturation transition introduces some uncertainty in the estimation of the accumulation thickness using the fp attribute. In other words, if the existence of a transition is not taken into account, for a given value of fp the thickness of the accumulation, as given by equation (8), is underestimated leading also to an incorrect estimation of its associated CO 2 volume.
A more extended and insightful study about the reflectivity and characterization of CO 2 saturation transitions was published by us in [13] and [14]. In those works, a comparative analysis of the results obtained with the layered saturation model and those obtained using the classic Wolf ramp [10], [11] was also presented.

Conclusions
In this chapter we presented a review on the theoretical seismic reflectivity function for carbon dioxide accumulations. For monitoring purposes, emphasis was placed on the study of its correlation with parameters of interest, such as accumulation thickness and overall saturation. This study is intended to predict quantitatively the changes that can be expected in seismic time lapse data analysis. Although for the examples we calibrated our rock physics models with information of the well known Sleipner field, similar analysis can be implemented for any other geological site, including the case of CO 2 -oil, or three phase fluid saturation with a proper adaptation of the physical models.
Given that the analysis of seismic AVA is a widely used tool in reservoir geophysics, using the classic Shuey's approximation, we modeled the behavior of intercept, gradient and curvature attributes under variable saturation and physical conditions. From our results we conclude that a monotonic increase in magnitude for the intercept can be expected in time lapse surveys, with significant variations with respect to the pre-injection state. This parameter shows clearly the decrease in the acoustic impedance in the sandstone for increasing CO 2 saturation.
The gradient is expected to decay strongly for low CO 2 saturations showing a different trend for saturations in the intermediate range, particularly for gaseous conditions. The curvature parameter is not expected to bring much more information about saturation, since in most cases its variations are smaller than those observed in the intercept and gradient. We have also shown that these AVA coefficients do not depend strongly on the choice of the equation of state to model the physical properties of CO 2 under in-situ reservoir conditions.
A loss of sensitivity on the AVA parameters, due to the stabilization of the seismic velocities, is characteristic of middle to high CO 2 saturations. This imposes a physical limitation on the monitoring of the CO 2 saturation. Our results suggest that a combined seismic time lapse analysis of the intercept and gradient AVA coefficients may be useful, not only to constrain the saturation state, but also to distinguish the physical state of the CO 2 accumulated below the caprock.
Taking into account that the vertical scale lengths in a layered medium provide physical constraints on the behavior of the propagating waves, we introduced this parameter in the computations of the reflectivity for a better characterization of this attribute. This is an issue of particular importance when the thicknesses of the CO 2 accumulations are small compared to the seismic wavelengths, a common situation in geologic repositories. The interference between multiple reflections gives rise to reflectivity dispersion effects and periodicities which can bring useful information. In this context, the use of modern spectral decomposition techniques plays an important role for the study of the amplitudes of reflected waves at different frequencies giving rise to amplitude versus frequency analysis, AVF. As shown in the examples, the periodicity of the reflectivity and the location of the first peak frequency can bring useful information about the thicknesses of the CO 2 bearing layers. In this sense, the utilization of the theoretical models presented here can be very useful for interpreting the observations and predicting peak frequencies. Our model estimates also show that the sensitivity of the reflectivity to saturation is frequency dependent, showing that low and high reflectivity bands can be expected. The high reflectivity zones show at the same time significant sensitivity to CO 2 saturation, a correlation that can be exploited for monitoring purposes.
We also analyzed the existing correlation between peak frequencies and thickness, to illustrate the utility of this spectral attribute in CO 2 accumulations. Our modeling results also indicate that the presence of saturation transitions below a main accumulation is a source of some inaccuracy in thickness and volume determination.
Regarding the potential utilization of these models and attributes (AVA coefficients, AVF, peak frequencies) using real seismic data to monitor changes in the CO 2 repository, some precautions must be taken to obtain meaningful estimations: • it is very important to use high quality time-lapse seismic data, with good repeatability and frequency content as wide as possible; • given that the amplitude spectrum of the seismic data is a combination of the wavelet spectrum and the reflectivity of the layers, it is very important for the processing sequence to include amplitude spectral balancing to remove wavelet overprint; • correct identification and isolation of a window containing the seismic reflections associated with the top of CO 2 accumulations is necessary; • the determination of peak frequencies and their time-lapse variations using an appropriate time-frequency decomposition on the selected seismic window is necessary.
In this way the theoretical results presented in this work can be tested for real case situations and may become useful tools for carbon dioxide monitoring problems.