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

__notoc__
{{OR/17/002}}
The seismic responses of different reservoirs to CO<sub>2</sub> 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<sub>2</sub> 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<sub>2</sub> 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<sub>2</sub> 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<sub>2</sub> 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<sub>2</sub> layers===
Interpretation in DiSECCS focussed on the topmost CO<sub>2</sub> 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<<ref name="Furre2015">>Furre, A K, Kiær, A, & Eiken, O. 2015. CO<sub>2</sub>-induced seismic time shifts at Sleipner, Interpretation, 3(3), 23–35.</ref>) 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.

[[Image:OR17002fig1.jpg|thumb|center|500px|'''Figure 1''' a) Map of the topmost CO<sub>2</sub> 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<sub>2</sub> 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<sub>2</sub> plus minor CH<sub>4</sub>) are uncertain, and, in particular, mixing scales (uniform or patchy) of the injected fluid and the aquifer brine are not well understood.

[[Image:OR17002fig2.jpg|thumb|center|500px|'''Figure 2''' Seismic velocity for the Utsira Sand as a function of CO<sub>2</sub> 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<sub>2</sub> layer that is essentially independent of the rock physics and its associated uncertainties. Key to this is understanding the geometry of the CO<sub>2</sub>-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<sub>2</sub> layer (Figure 3).

[[Image:OR17002fig3.jpg|thumb|center|500px| '''Figure 3''' a) Relationship between top layer topography (ΔH<sub>E</sub> ,  ΔT<sub>E</sub> ) and time-shifts  of the base layer reflection (velocity pushdown ΔT<sub>P</sub>) assuming a horizontal CO<sub>2</sub> – water contact and no wavelet interference effects. b) Seismic section through a ridge showing schematic measurement of ΔT<sub>E</sub> and ΔT<sub>P</sub> 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:

{|
| [[Image:OR17002equation1.jpg|frameless|100px]]
| Equation 1.
|}

Where:<br>
ΔT<sub>E</sub> = time-shift at top layer reflection<br>
ΔT<sub>P</sub> = time-shift at base layer reflection<br>
V<sub>O</sub> = overburden velocity<br>
V<sub>R</sub> = CO<sub>2</sub> layer velocity

Rearranging Equation 1 gives an expression for V<sub>R</sub> as a function of the overburden velocity and the time-shifts at the top and base of the CO<sub>2</sub> layer:

{|
| [[Image:OR17002equation2.jpg|frameless|100px]]
| Equation 2.
|}

ΔT<sub>E</sub> and ΔT<sub>p</sub> can readily be measured from seismic data and V<sub>o</sub> can be estimated from borehole and seismic processing velocities, so V<sub>R</sub> can be calculated. Early results (Figure 4) were broadly in line with the rock physics and indicative of moderate to high CO<sub>2</sub> saturations (Chadwick et al. 2016<ref>Chadwick, R A, Williams, G A, & White, J C. 2016. High resolution imaging and characterisation of a CO<sub>2</sub> layer at the Sleipner CO<sub>2</sub> Storage operation using time-lapse seismics. First Break, 34, 79–87.</ref>).

[[Image:OR17002fig4.jpg|thumb|center|600px|'''Figure 4''' Calculated layer velocity distributions measured on multiple seismic traces for different ridge elevation (''Δ''T<sub>E</sub>) thresholds a) ''Δ''T<sub>E</sub> > 2 ms b) ''Δ''T<sub>E</sub> > 3 ms c) ''Δ''T<sub>E</sub> > 4&nbsp;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<sub>2</sub>-water contact level and also a means of correcting for interference induced time-shifts in the reflected wavelets.

[[Image:OR17002fig5.jpg|thumb|center|400px|'''Figure 5''' Layer velocities calculated for the topmost CO<sub>2</sub> 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<sub>2</sub> 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<ref name="Williams2013">Williams, G A, & Chadwick, R A. 2013. Quantitative seismic analysis of a thin layer on CO<sub>2</sub> in the Sleipner injection plume. Geophysics. 77(6), 245–256.</ref>), 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).

[[Image:OR17002fig6.jpg|thumb|center|400px|'''Figure 6''' Left: Temporal spacings extracted from the central part of the layer shown in (a), where the directly measured temporal spacing exceeds 8&nbsp;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).

[[Image:OR17002fig7.jpg|thumb|center|400px|'''Figure 7''' Temporal thicknesses (milliseconds) of the CO<sub>2</sub> 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<sub>2</sub> 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)<ref name="Lyle2005">Lyle, S, Huppert, H E, Hallworth, M, Bickle, M, & Chadwick, A. 2005. Axisymmetric gravity currents in a porous medium. ''Journal of Fluid Mechanics'' 543, 293–302.</ref>. Numerical models included full 3-D and high resolution 2D axisymmetric geometry.

[[Image:OR17002fig8.jpg|thumb|center|400px|'''Figure 8''' 2D axisymmetric reservoir model showing the effect of injecting CO<sub>2</sub> at 48°C into a cooler reservoir. The reservoir temperature at the injection point (labelled) is ~35°C. (a) Simulation mesh, (b) CO<sub>2</sub> saturation in the reservoir, (c) thermal anomaly (°C), and (d) CO<sub>2</sub> density (kgm<sup>-3</sup>). The simulation time-step corresponds to the 2006 seismic monitor survey. The radial equivalent extents of each seismically observed CO<sub>2</sub> layer in the plume are marked with a white square in (b).]]

The spread of the CO<sub>2</sub> 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<sub>2</sub> 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<sub>2</sub> 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<sub>2</sub> mobility, at least in the axial part of the plume (Figure 8c, d).

By including higher reservoir permeability and higher CO<sub>2</sub> temperature we concluded (Williams & Chadwick 2017<ref>Williams, G, & Chadwick R A. 2017. An improved history-match for layer spreading within the Sleipner plume including thermal propagation effects. Energy Procedia.

</ref>) 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<sub>2</sub>.

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<sub>2</sub> injection scenarios and also comparing them with analytical models of a thin spreading CO<sub>2</sub> layer (Williams et al. 2018<ref>Williams, G A, Chadwick, R A, & Vosper, H. 2018. Some thoughts on Darcy-type flow simulation for modelling underground CO<sub>2</sub> storage based on the Sleipner CO<sub>2</sub> storage operation. ''International Journal of Greenhouse Gas Control'', 68 (2018), 164–175.</ref>). 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.

[[Image:OR17002fig9.jpg|thumb|center|200px|'''Figure 9''' Effect of injecting warm CO<sub>2</sub> into the topmost sand body. (a) Ambient temperature distribution at the top of the reservoir. (b) CO<sub>2</sub> distribution at the time of the 2010 seismic survey, assuming the CO<sub>2</sub> enters the top sand body at the background reservoir temperature. (c) Thermal anomaly caused by CO<sub>2</sub> entering the top sand body at a temperature of 37°C. (d) Resulting CO<sub>2</sub> distribution at the time of the 2010 seismic monitor survey. Black polygon denotes extent of the CO<sub>2</sub> measured on the 2010 survey. White circle denotes feeder chimney through which the CO<sub>2</sub> enters the topmost sand from below. Temperature and thermal anomaly in °C, plume distribution shown as CO<sub>2</sub> saturation.]]

===Saturation&nbsp;—&nbsp;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<sub>2</sub> plume and more widespread pressure changes in the whole reservoir. Thus high frequency tuning occurs around the thin layers of CO<sub>2</sub> 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<sub>2</sub>, ponding and spreading radially beneath an impermeable topseal. Analytical modelling (Lyle et al. 2005<ref name="Lyle2005"></ref>) 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<sub>2</sub> 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<sub>2</sub> 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<sub>2</sub> saturation in a thin layer is a generic and critical problem for CO<sub>2</sub> monitoring (e.g. Arts et al., 2004<ref>Arts, R, Eiken, O, Chadwick, A, Zweigel, P, van der Meer, B, & Kirby, G. (2004). Seismic monitoring at the Sleipner underground CO<sub>2</sub> storage site (North Sea). Geological Society, London, Special Publications, 233(1), 181–191.</ref>; Chadwick et al. 2004<ref>Chadwick, R A, Arts, R, Eiken, O, Kirby, G A, Lindeberg, E, & Zweigel, P. 2004. 4D seismic imaging of an injected CO<sub>2</sub> bubble at the Sleipner Field, central North Sea. In: Davies, R J, Cartwright, J A, Stewart, S A., Lappin, M, & Underhill, J R. (eds) 3D Seismic Technology: Application to the Exploration of Sedimentary Basins. Geological Society, London, Memoir 29, 305–314, The Geological Society of London 2004.</ref>). 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<ref>Müller, T M, Gurevich, B, & Lebedev, M. 2010. Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks—A review. Geophysics, ''75''(5), 75A147–75A164.</ref>; Chapman et al., 2002<ref name="Chapman2002">Chapman, M, Zatsepin, S V. & Crampin, S. 2002. Derivation of a microstructural poroelastic model. ''Geophysical Journal International'', 151(2), 427–451.</ref>). 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''.<br>
: b) ''Fast and efficient wavefield modelling for the effects of dispersion in complex velocity models''.<br>
: c) ''Numerical simulation of the effect of dispersion on saturation estimation methods''.<br>
: 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)<ref name="Chapman2002"></ref>, 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<sub>2</sub> storage in Snøhvit (Grude et al., 2013<ref>Grude, S, Landrø, M, & Osdal, B. 2013. Time-lapse pressure–saturation discrimination for CO<sub>2</sub> storage at the Snøhvit field. ''International Journal of Greenhouse Gas Control'', 19, 369–378.</ref>), Sleipner (Carcione et al., 2006<ref>Carcione, J M, Picotti, S, Gei, D, & Rossi, G. 2006. Physics and seismic modeling for monitoring CO<sub>2</sub> storage. Pure and Applied Geophysics, 163(1), 175–207.</ref>) and Nagaoka in Japan (Lumley, 2010<ref>Lumley, D. 2010. 4D seismic monitoring of CO<sub>2</sub> sequestration. ''The Leading Edge'', 29(2), 150–155.</ref>).

Good agreement is seen between the theory and Brie’s model (Papageorgiou and Chapman, 2015<ref>Papageorgiou, G, & Chapman M. 2015. Multi-fluid squirt flow and Hysteresis Effects on the Bulk Modulus-Water Saturation relationship. ''Geophysical Journal International''. 203(2), 814–817. [http://gji.oxfordjournals.org/content/203/2/814.short http://gji.oxfordjournals.org/content/203/2/814.short]</ref>), 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.

[[Image:OR17002fig10.jpg|thumb|center|400px|'''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<sub>2</sub> (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).

[[Image:OR17002fig11.jpg|thumb|center|300px|'''Figure 11''' Seismic velocity (Vp) variation with supercritical CO<sub>2</sub>-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).]]

[[Image:OR17002fig12.jpg|thumb|center|300px|'''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<ref>Jin, Z, Chapman, M, Wu, X, & Papageorgiou, G. 2016. Estimating gas saturation in a thin layer by using frequency‐dependent amplitude versus offset modelling. Geophysical Prospecting.</ref>). 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====
<u>Matching pursuit methods</u><br>
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<sub>2</sub> 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<sub>2</sub> -saturated sand layer.

<u>Elastic dictionary</u><br>
In this approach we build the dictionary using a wedge model consisting of a CO<sub>2</sub> -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)<ref name="Furre2015"></ref>. The dictionary is built using elastic-only reflectivity on the shale-CO<sub>2</sub> and CO<sub>2</sub>-brine interfaces and different fluid mixing laws are considered for the calculation of the velocities.

The different fluid averaging laws&nbsp;—&nbsp;Wood, Brie and Voigt&nbsp;—&nbsp;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)<ref name="Chadwick2006">Chadwick, A, Arts, R, Eiken, O, Williamson, P, & Williams, G. 2006. Geophysical monitoring of the CO2 plume at Sleipner, North Sea. In Advances in the geological storage of carbon dioxide (pp. 303–314). Springer Netherlands.</ref>. 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)<ref name="Chadwick2006"></ref> and one by Williams & Chadwick (2013)<ref name="Chadwick2013">Williams, G A, & Chadwick, R A. 2013. Quantitative seismic analysis of a thin layer on CO<sub>2</sub> in the Sleipner injection plume. ''Geophysics''. 77(6), 245–256.</ref>. The importance of our method lies in that it is blind so it provides a third, independent, way to distinguish CO<sub>2</sub> layer thickness.

<u>Frequency-dependent dictionary</u><br>
In this approach the dictionary is built on the assumption the reflectivity of the sand is dispersive which means that both the shale-CO<sub>2</sub> and CO<sub>2</sub>-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)<ref name="Furre2015"></ref>, 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&nbsp;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.

[[Image:OR17002fig13.jpg|thumb|center|400px|'''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)<ref name="Chadwick2006"></ref> and the dashed black line shows the estimate using the spectral method of Williams and Chadwick (2013)<ref name="Williams2013"></ref>.]]

[[Image:OR17002fig14.jpg|thumb|center|400px|'''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.]]

[[Image:OR17002fig15.jpg|thumb|center|400px|'''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.]]

==References==
{{reflist}}

[[Category: OR/17/002 DiSECCS - Diagnostic Seismic Toolbox for Efficient Control of CO2 Storage: Final Summary Report: Work Packages 1-4 | 05]]