OR/17/002 Work Package 2: Monitoring to understand detailed flow processes and improve in situ quantification

Chadwick, R A. 2017. DiSECCS Final Summary Report. Work Packages 1 – 4 British Geological Survey Internal Report, OR/17/002.

The seismic responses of different reservoirs to CO₂ injection involve a number of process‐related trade-offs. High quality (thick, permeable and mechanically compliant) reservoirs tend to have a large seismic sensitivity and relatively high potential for the spatial and temporal resolution of thin, spreading CO₂ layers and the characterisation of layer properties. But they show rather small pressure increases, close to or beneath seismic detection limits. On the other hand, lower quality reservoirs tend to have a decreased seismic sensitivity with less potential for CO₂ layer characterisation but, conversely, with larger and more seismically tractable pressure increases.

DiSECCS had access to monitoring datasets which covered a range of reservoir quality: notably from the Utsira Sand at Sleipner, a thick homogeneous reservoir; the Stø Formation at Snøhvit also quite homogeneous; and the much more structurally and stratigraphically complex Tubåen reservoir at Snøhvit. Because of this we were able to address a number of different issues related to understanding injection processes, including characterisation of thin CO₂ layers, discrimination between fluid saturation and pressure effects and modelling layer flow.

Integrated very high resolution flow modelling and seismic analysis (Task 2.1)

Work in Task 2.1 focussed on very accurate measurement of CO₂ layer morphology and temporal thickness, integrated with geological constraints and fluid flow modelling. The high-resolution 2010 seismic survey from Sleipner comprised the key dataset for this analysis, with the 2012 seismic data from Snøhvit also providing useful insights.

Morphology of thin CO₂ layers

Interpretation in DiSECCS focussed on the topmost CO₂ layer in the injected Sleipner plume which has accumulated beneath the reservoir topseal relief. The high-resolution 2010 seismic data have, in places, achieved explicit 3D imaging of the upper and lower reflective surfaces of the layer (Figure 1). It has been noted previously (Furre et al. 2015<) that interference between the top and base layer reflections produces subtle time-shifts which can be strongly diagnostic of layer temporal thickness, and this approach to forensic interpretation of layer reflectivity and travel-times has been taken forward in DiSECCS.

Figure 1 a) Map of the topmost CO₂ layer in 2010 showing temporal thicknesses and line of section b) section through topmost layer showing explicit resolution (separation) of top and base layer reflections (blue and orange respectively) beneath two north-trending ridges at the upper surface of the reservoir.

Layer velocity

CO₂ layer velocities at Sleipner can be estimated from rock physics (Figure 2), but they remain rather poorly-constrained. This is because only a single core is available which might not be representative of the reservoir as a whole, properties of the injected fluid (CO₂ plus minor CH₄) are uncertain, and, in particular, mixing scales (uniform or patchy) of the injected fluid and the aquifer brine are not well understood.

Figure 2 Seismic velocity for the Utsira Sand as a function of CO₂ saturation, for uniform (Reuss average) and patchy (Voigt average) fluid mixing scales.

Small time-shifts were used to obtain an empirical determination of the seismic velocity of the topmost CO₂ layer that is essentially independent of the rock physics and its associated uncertainties. Key to this is understanding the geometry of the CO₂-water contact. This forms the base layer reflection which is subject to small time-shifts depending on the topography of the overlying topseal and the thickness and velocity of the CO₂ layer (Figure 3).

Figure 3 a) Relationship between top layer topography (ΔH_E , ΔT_E ) and time-shifts of the base layer reflection (velocity pushdown ΔT_P) assuming a horizontal CO₂ – water contact and no wavelet interference effects. b) Seismic section through a ridge showing schematic measurement of ΔT_E and ΔT_P on five notional seismic traces (dashed lines).

With a perfectly resolved seismic image (i.e. for a seismic ‘spike’ with no interference between the top and base layer wavelets) there is a simple relationship between time-shifts due to topographic variation at the top layer reflection and time-shifts due to velocity pushdown at the base layer reflection:

Equation 1.

Where:
ΔT_E = time-shift at top layer reflection
ΔT_P = time-shift at base layer reflection
V_O = overburden velocity
V_R = CO₂ layer velocity

Rearranging Equation 1 gives an expression for V_R as a function of the overburden velocity and the time-shifts at the top and base of the CO₂ layer:

Equation 2.

ΔT_E and ΔT_p can readily be measured from seismic data and V_o can be estimated from borehole and seismic processing velocities, so V_R can be calculated. Early results (Figure 4) were broadly in line with the rock physics and indicative of moderate to high CO₂ saturations (Chadwick et al. 2016).

Figure 4 Calculated layer velocity distributions measured on multiple seismic traces for different ridge elevation (ΔT_E) thresholds a) ΔT_E > 2 ms b) ΔT_E > 3 ms c) ΔT_E > 4 ms

Subsequently we developed a number of top layer topographic models, with a range of layer velocities, and computed synthetic seismic datasets to cover all likely topseal topography scenarios. From these we developed a methodology of identifying a more accurate CO₂-water contact level and also a means of correcting for interference induced time-shifts in the reflected wavelets.

Figure 5 Layer velocities calculated for the topmost CO₂ layer at Sleipner in 2010 showing uncorrected velocity values and values corrected for top/base reflection interference effects.

Revised velocity estimates incorporating the synthetic model constraints (Figure 5) are rather lower than the earlier values, with a slightly smaller scatter.

Temporal layer thicknesses

In parallel with the above, a comparative study of multiple approaches to obtaining temporal thicknesses of the topmost CO₂ layer at Sleipner was carried out. These included direct measurement of the resolved layer from high resolution data (as above), frequency analysis using the smoothed pseudo Wigner-Ville distribution (SPWVD) spectral decomposition tool (Williams & Chadwick 2013), modified amplitude analysis, structural analysis of topseal relief and time-shifts both of the layer reflections (edge-pull-up) and of reflections beneath the layer (velocity push-down) (Figure 6).

Figure 6 Left: Temporal spacings extracted from the central part of the layer shown in (a), where the directly measured temporal spacing exceeds 8 ms on the 2010IP data. Middle: Cross-plots of temporal spacing derived from each analysis technique against corrected measured temporal spacing. Right: Normalised histograms of the difference between the spacing derived at each seismic trace using the particular technique and the corrected measured temporal spacing. The green lines indicate an idealised 1:1 correlation whilst the red line shows the actual linear best fit with an intercept through the origin.

Results indicate that, depending on layer thickness, different methods give the best estimates of temporal thickness and an optimal strategy for layer volume estimation could utilise a number of methods in combination. For example measured temporal spacings would be used in the central part of the layer, above the tuning thickness, spectral decomposition would be used down to the effective resolution limit and reflection amplitudes for the thinnest parts of the layer out to the layer edge. A quantification of the topmost layer at Sleipner, integrating a number of techniques, has been carried out (White et al. submitted).

Figure 7 Temporal thicknesses (milliseconds) of the CO₂ plume in the Stø reservoir at Snøhvit in 2011 and 2012 a) from direct measurement and amplitude scaling b) from spectral decomposition.

A similar approach was taken to analyse the 2012 seismic dataset from Snøhvit which images second phase of CO₂ injection, into the Stø Formation (White et al. submitted). Here temporal layer spacings measured directly in the (thicker) central part of the plume, combined with temporal spacings derived from the SPWVD tuning frequencies and amplitude analysis, allows reliable mapping of temporal layer thickness across the plume extent (Figure 7).

Flow models

Fluid flow modelling work was carried out in conjunction with the seismic analysis. A number of numerical flow simulators were utilised, including TOUGH2, PFLOTRAN, Eclipse 100 and Eclipse 300, and analytical models of thin layer spreading were also developed using the axisymmetric gravity current formulation of Lyle et al. (2005). Numerical models included full 3-D and high resolution 2D axisymmetric geometry.

Figure 8 2D axisymmetric reservoir model showing the effect of injecting CO₂ at 48°C into a cooler reservoir. The reservoir temperature at the injection point (labelled) is ~35°C. (a) Simulation mesh, (b) CO₂ saturation in the reservoir, (c) thermal anomaly (°C), and (d) CO₂ density (kgm^-3). The simulation time-step corresponds to the 2006 seismic monitor survey. The radial equivalent extents of each seismically observed CO₂ layer in the plume are marked with a white square in (b).

The spread of the CO₂ plume at Sleipner has been a topic of much discussion in recent years, and a better understanding of this was a key research objective. Our focus was on the accuracy of the reservoir geological characterisation and on the thermal setting. Firstly, by properly integrating core measurements in the Utsira Sand with the geophysical log data, we obtained an improved assessment of reservoir permeability. Secondly, Statoil have estimated that the CO₂ is injected at about 13°C above reservoir temperature, but this situation had not previously been modelled. Our simulation of the thermal effects of injecting ‘warm’ CO₂ at above the ambient reservoir temperature (Figure 8) suggests upward propagation of the thermal anomaly to the reservoir top with a significant effect on CO₂ mobility, at least in the axial part of the plume (Figure 8c, d).

By including higher reservoir permeability and higher CO₂ temperature we concluded (Williams & Chadwick 2017) that numerical simulations based on the assumption of Darcy flow can satisfactorily account for the observed rapid layer spreading (Figure 9). In fact the principle factor is the revised geological reservoir parameters, rather than the warmer CO₂.

In parallel with this work, we have carried out a very accurate flow simulation code comparison, comparing the relative performance of a number of numerical flow simulators on a range of carefully matched CO₂ injection scenarios and also comparing them with analytical models of a thin spreading CO₂ layer (Williams et al. 2018). A key conclusion is that the simulators all produced very similar results, with differences between the different codes less than interpretive uncertainty in the monitoring data.

Figure 9 Effect of injecting warm CO₂ into the topmost sand body. (a) Ambient temperature distribution at the top of the reservoir. (b) CO₂ distribution at the time of the 2010 seismic survey, assuming the CO₂ enters the top sand body at the background reservoir temperature. (c) Thermal anomaly caused by CO₂ entering the top sand body at a temperature of 37°C. (d) Resulting CO₂ distribution at the time of the 2010 seismic monitor survey. Black polygon denotes extent of the CO₂ measured on the 2010 survey. White circle denotes feeder chimney through which the CO₂ enters the topmost sand from below. Temperature and thermal anomaly in °C, plume distribution shown as CO₂ saturation.

=Saturation — pressure discrimination by frequency analysis of thin and thick layers

Working on the Snøhvit time-lapse data, White et al. (2015) showed that for injection into the Tubåen reservoir, frequency analysis could be used to discriminate between fluid saturation changes associated with the CO₂ plume and more widespread pressure changes in the whole reservoir. Thus high frequency tuning occurs around the thin layers of CO₂ in the plume, contrasted with much lower frequency tuning associated with fluid pressure changes propagating across the full reservoir thickness.

Similar analysis of the more recent injection into the Stø Formation (White et al. in submission), shows no evidence of low frequency tuning away from the injection plume, which suggests a lack of significant pressure propagation into the reservoir. The seismic anomaly is more simple, and consistent with a cone-shaped plume formed by buoyancy-driven upward advection of CO₂, ponding and spreading radially beneath an impermeable topseal. Analytical modelling (Lyle et al. 2005) of radially-symmetric gravity currents shows that layer thickness varies with rate of injection, so the observed constancy is consistent with the roughly uniform rate of injection.

Novel methods for layer characterisation (Task 2.2)

This work package has addressed the seismic characterization of CO₂ saturation in thin layers. We have developed new theory for calculating frequency-dependent velocity and attenuation in rocks saturated by multiple fluids. The key advance is to describe the 'squirt' and 'patch' effects in a common theoretical framework. The new model reproduces existing work in suitable limits and demonstrates that patch effects can be produced by complex pore-scale fluid distributions, so that interpretation of data in terms of explicit patch sizes is problematic. The model has been calibrated against experimental data obtained in WP3, and provides a compelling explanation of velocity and attenuation as a function of CO₂ saturation. Combination of our rock physics models with a matching pursuit method allows application to layer characterization at Sleipner. Using only the frequency-independent theories the method produces results which can be compared to previous estimates using independent techniques. Repeating the analysis with the new modelling, we demonstrate the potential importance of including dispersion and attenuation for saturation estimation.

Estimation of CO₂ saturation in a thin layer is a generic and critical problem for CO₂ monitoring (e.g. Arts et al., 2004; Chadwick et al. 2004). Approaching this problem requires knowledge of the background velocity model for the reservoir and surrounding rocks, together with an understanding of the rock physics which link velocities and density to saturation. Finally, we need to be able to perform wavefield modelling which takes account of the amplitude versus angle response, travel-times through the saturated layers and the interference between top and bottom reflections.

The rock physics relevant to wave propagation through rocks saturated with multiple fluids is not fully understood. We expect multi-fluid saturation to give rise to strong dispersion and attenuation (e.g. Müller et al, 2010; Chapman et al., 2002). Such effects would have potentially important implications for estimation of fluid saturation within a thin layer.

Being able to address the issues of dispersion and attenuation in the context of practical monitoring of CCS relies on solution to a number of problems:

a) Reliable modelling tools for the calculation of frequency-dependent velocity and attenuation.

b) Fast and efficient wavefield modelling for the effects of dispersion in complex velocity models.

c) Numerical simulation of the effect of dispersion on saturation estimation methods.

d) Demonstration of the techniques on suitable field datasets

We have addressed each point within the DiSECCS research programme, and treat them in turn below.

Reliable modelling of frequency-dependent velocity and attenuation

The rock physics relevant to wave propagation through rocks saturated with multiple fluids is not fully understood. Modelling generally makes use of the patchy or uniform saturation theories, which are highly idealized and neglect potentially important effects. The difference in behaviour between the uniform and patch theories stems from a difference in physical assumptions-under the uniform model a single fluid pressure always obtains in each fluid and at each point in the pore space, while under the patch model there can be two such pressures, one for each fluid.

The patch effect itself can give rise to frequency-dependent velocity and attenuation, but such effects can often be understood well within the context of the 'squirt-flow' theory. Combining the squirt and patch theories is a key research objective, which we have addressed by developing a new theory which combines the 'squirt' and 'patch' concepts.

The basis of the work is to extend previous analysis by assuming that induced pore-pressures are systematically different between the two saturating fluids. This can be due either to membrane effects or differences in fluid distributions. The difference between the fluid pressures is captured by a single new capillary pressure parameter, and with the addition of this relation we are able to re-derive the earlier squirt flow models based on Chapman et al. (2002), obtaining a new theory which contains both squirt and patch effects under a consistent parameterization.

The new theory reproduces key theoretical approaches, such as Gassmann-Wood and single fluid squirt flow theory in appropriate limits and predicts that the characteristic frequency where squirt flow effects are observed depends on relative permeability. A particularly interesting limit is the low frequency limit in which the behaviour is dependent only on the capillary pressure parameter. These results can be compared directly with the predictions of the so-called 'Brie model', which is widely used in fluid substitution problems related to CO₂ storage in Snøhvit (Grude et al., 2013), Sleipner (Carcione et al., 2006) and Nagaoka in Japan (Lumley, 2010).

Good agreement is seen between the theory and Brie’s model (Papageorgiou and Chapman, 2015), so that we have in effect given a theoretical derivation of what was an empirical model. This builds confidence in our approach, as well as allowing us to enhance existing interpretations based on the Brie model. A schematic representation of the various limits of the new theory is shown in Figure 10.

Figure 10 The partial-saturation squirt-flow model (upper left) developed here incorporates the effects of patchy saturation by letting the capillary pressure parameter be less than 1. It also incorporates the frequency-dependent effects of the squirt flow model of Chapman et al. (2002)ref name="Chapman2002"></ref>. Its limits depend on which of these parameters (frequency, capillary pressure parameter) is switched off.

We have also worked with WP3 to calibrate our new model against ultrasonic velocity and attenuation experimental data obtained at different partial saturation of supercritical CO₂ (Figures 11 and 12) The simultaneous fit of attenuation and ultrasonic velocity versus saturation means we constrained rock physics parameters determining the magnitude of squirt flow (crack density, characteristic frequency).

Figure 11 Seismic velocity (Vp) variation with supercritical CO₂-brine saturation and the fit of the squirt flow model for the ultrasonic velocity measurements of the Utsira Sand. The solid line is Gassmann’s model with a Wood-averaged fluid, the dot-dashed line with a Voigt–averaged fluid. The dashed line is the best fit of our model to the data and the dotted line is its low frequency limit (data from WP3).

Figure 12 Attenuation variation with water saturation for the Utsira Sand. The fitted line is the attenuation corresponding to the squirt flow (dashed line) model shown in Figure 11 (data from WP3).

Convolutional based modelling incorporating the effects of frequency-dependent reflectivity

Realistic application of the models described above to the characterization of thin layers requires some form of seismic inversion. Most commonly employed seismic inversion is based on the convolutional model of a seismic trace. The convolutional model is computationally very efficient and has the key advantage that it simulates ideal data after processing. More sophisticated methods such as finite difference and reflectivity provide more rigorous solutions to the wave equation, but they model shot gathers and this output must itself be processed before it can be compared to seismic data. This requirement impedes the application of such methods to inversion of stacked data or pre-stack gathers.

Methods to incorporate frequency-dependent velocities and attenuations into finite difference and reflectivity methods are well established, but convolution-style modelling incorporating the effects of frequency-dependent reflectivity has been less well studied. To address this issue, we have developed an extended convolutional model which can compute synthetic seismograms in the time-angle domain on the basis of our developed rock physics models (see Jin et al., 2016). Unlike the conventional convolutional model, which accepts Vp, Vs and density logs as inputs, we also require porosity and fluid saturation to be specified, which opens the way to potentially inverting for saturation as a function of depth.

In the simplest case of a thin layer under normal incidence, the seismic response is the result of a complex interplay between thin-layer tuning, frequency-dependent reflectivity and phase effects on the wavelet resulting from the effect of attenuation on the reflection coefficient.

Application of the modelling to field data

Matching pursuit methods
Our approach to the analysis of the Sleipner dataset is based on combining matching pursuit methods and rock physics to detect thin layers and characterize saturation and thickness. Matching pursuit methods are a common means of decomposing seismic data into a sum of user-defined basis functions, specified in a so-called 'dictionary'. For our analysis, we specify a dictionary containing forward-modelled responses for various layer thicknesses, based on our rock physics model, so that the matching pursuit algorithm searches for these responses in the data. In this way, the algorithm both detects the presence of the thin layer, and finds the best fitting combination of thickness and saturation.

We applied this algorithm with the aim to determine the thickness of the layer of CO₂ visible in the Sleipner 2010 time-lapse data within inlines 1871–1877 and crosslines 1200–1220. We can construct our dictionary choosing only elastic rock physics models, or we can include the frequency-dependent rock physics models developed during DiSECCS for the reflectivity of the CO₂ -saturated sand layer.

Elastic dictionary
In this approach we build the dictionary using a wedge model consisting of a CO₂ -saturated sand of varying thickness encased between a shale caprock and a brine-saturated sand. Suitable values for the velocity and density were taken from Furre et al. (2015). The dictionary is built using elastic-only reflectivity on the shale-CO₂ and CO₂-brine interfaces and different fluid mixing laws are considered for the calculation of the velocities.

The different fluid averaging laws — Wood, Brie and Voigt — are substituted into Gassman’s model and we then build the dictionary for each fluid mixing scenario based on the wedge model. By fixing the saturation, we invert for layer thickness. The method produces results comparable to those obtained using the overburden method introduced by Chadwick et al. (2006). In Figure 13 we show the result of this inversion together with published estimates of layer thickness in that location: one by Chadwick et al. (2006) and one by Williams & Chadwick (2013). The importance of our method lies in that it is blind so it provides a third, independent, way to distinguish CO₂ layer thickness.

Frequency-dependent dictionary
In this approach the dictionary is built on the assumption the reflectivity of the sand is dispersive which means that both the shale-CO₂ and CO₂-brine reflectivities are frequency dependent. Consequently, at each interface there are phase effects dependent on the attenuation of the layer encoded in the dictionary.

We calculated synthetic seismograms for a variety of layer thicknesses and saturations both with and without frequency-dependent reflectivity, for a Ricker wavelet with a central frequency of 35 Hz. Measuring the amplitude versus the time thickness for various saturations gave rise to the tuning curves shown in Figure 14. The curves for the elastic case are given on the left, and these are consistent with the work of Furre et al. (2015), but we note that the behaviour changes for the frequency-dependent case, due to the phase effects on the reflection coefficient.

This motivates us to compare inversions for saturation and thickness using both dispersive and non-dispersive models. We found that the phase effects have a modest impact on the inverted layer thicknesses (generally less than 2 m difference), but there is a more pronounced impact on the saturation estimation. Figure 15 shows the residuals of the elastic and frequency dependent inversion. The minimum for the elastic and dispersive dictionaries are different suggesting that it is important to consider dispersive effects for the estimation of saturation in a thin layer.

Figure 13: Estimates of layer thickness for inlines 1871 to 1877 using matching pursuit under different fluid substitution models: Wood (red), Brie (green) and Voigt (blue). The solid black line denotes the estimate of thickness using the overburden method of Chadwick et al. (2006) and the dashed black line shows the estimate using the spectral method of Williams and Chadwick (2013).

Figure 14: Amplitude tuning curves for a Ricker wavelet of central frequency 35 Hz. The elastic response on the left is derived using Gassmann’s model with a Brie mixing law and on the right the dispersive rock physics theory has been used. The two methods give different tuning amplitudes and minima for different degrees of water saturation.

Figure 15 Sum of residuals of the matching pursuit inversion for the whole saturation range in 10% saturation steps. The residuals refer to the sum over the entire 7x21 cube of traces in the location of interest. The dispersive and elastic dictionaries have minima (red circles) at different saturations.