CO 2 -brine ﬂ ow-through on an Utsira Sand core sample: Experimental and modelling. Implications for the Sleipner storage ﬁ eld

Sleipner (North Sea) is the world ’ s ﬁ rst commercial-scale carbon capture and storage (CCS) project, active since 1996, with ∼ 17 million tonnes of CO 2 stored. The main reservoir, Utsira Sand, constitutes an ideal host formation of exceptionally high porosity-permeability and large lateral extent. However, the extensive seismic time-lapse, gravity and electromagnetic monitoring surveys deployed at Sleipner have not been well-supported by laboratory measurements. Here, we investigate the geophysical and geomechanical response of an Utsira core sample for the ﬁ rst time, using controlled in ﬂ ation/depletion cycles at variable CO 2 -to-brine fractional ﬂ ow rates. Ultrasonic P-wave velocities and attenuations are measured together with electrical resistivity (converted into CO 2 -saturation), along with continuous axial and radial strain monitoring. Ultrasonic velocity and at- tenuation data were simultaneously inverted and results extrapolated to ﬁ eld-scale seismic-frequencies using a new rock physics theory, which combines patchy ﬂ uid distribution and squirt ﬂ ow e ﬀ ects. It provides a velocity- saturation relationship of practical importance to CO 2 plume monitoring. Furthermore, by combining ultrasonic and deformation data, we report empirical relations between pore pressure changes and geomechanical e ﬀ ects in the reservoir, for di ﬀ erent saturation ranges. Our dataset complements and constrains existing geophysical monitoring surveys at Sleipner and, more generally, improves the understanding of shallow weakly-cemented sand reservoirs. Utsira Sand core sample during a CO 2 brine ﬂ ooding experiment. We measure longitudinal (primary) acoustic wave velocity and attenuation, and electrical resistivity together with a continuous record of axial and radial strains. The test covers seven drainage brine-CO 2 fractional ﬂ ow episodes simulating the transitional stages of increasing CO 2 saturation in a storage reservoir, plus an additional forced imbibition (brine ﬂ ow), after cessation of CO 2 injection. to investigate


Introduction
The Sleipner project has demonstrated the potential of carbon capture and storage (CCS) to be a realistic large-scale greenhouse gas mitigation technique. In 1996, Statoil and its Sleipner partners pioneered CO 2 storage at commercial-scale in Sleipner West field in the Norwegian North Sea (Baklid et al., 1996). Since then, around 1 Mt per year of CO 2 has been injected into the Utsira Sand, at a depth of ∼1020 m below sea level through a deviated well (Arts et al., 2004). There is no evidence of leakage above the host formation (Eiken et al., 2011).
The Utsira Sand is a regional saline aquifer comprising weakly-cemented sands of late Cenozoic age, overlain by a ∼700 m thick dominantly argillaceous overburden . The aquifer forms a high porosity (35-42%) high permeability (> 1D) low structural relief geological system Williams and Chadwick, 2012), some 200 m thick in the Sleipner area and ideal for CO 2 storage. Such a suitable combination of porosity-permeability and lateral extent (∼26.000 km 2 ) has resulted in very little pore pressure increase (< 0.1 MPa) over the last 16 years of continuous injection .
Comprehensive geophysical monitoring has been carried out on Sleipner: a set of eight 3D seismic surveys to 2014(Furre et al., 2015, three time-lapse gravity surveys and a trial of a Controlled Source Electromagnetic (CSEM) survey in 2008 (Eiken et al., 2011). The Pwave velocity (V p ) variation has successfully provided valuable information about subsurface CO 2 location, due to the velocity reduction of seismic waves travelling through the CO 2 -saturated layers (Chadwick et al., 2010). Quantification in situ of the CO 2 stored is more challenging and requires establishing an accurate relationship between V p and CO 2 saturation, which depends on fluid distribution patterns within the twophase fluid flow system (Yamabe et al., 2016).
Electromagnetic surveys represent an alternative tool for interpreting the movement of a CO 2 plume in saline aquifers due to strong T resistivity contrasts between brine and CO 2 (Alemu et al., 2013;Falcon-Suarez et al., 2016;Wang et al., 2009), and can potentially be used to quantify partial CO 2 saturations (Carrigan et al., 2013;Falcon-Suarez et al., 2017;Falcon-Suarez et al., 2016;Nakatsuka et al., 2010). However, they are of much lower resolution than the 3D seismic, with interpretation of the CO 2 plume at Sleipner further hampered by the lack of a baseline (pre-injection) survey (Park et al., 2013). Efforts are being made towards extracting further information from CSEM by combining seismic and electrical data, and also from rock physics modelling (Park et al., 2017).
Rock physics modelling makes use of observed velocity changes in seismic time-lapse data (e.g., Chadwick et al., 2010) by employing a rock physics based velocity-CO 2 saturation relationship to map the extent and saturation of the CO 2 plume. But the method carries a degree of uncertainty as the distribution of pore fluids affects the relationship. In the case of 'patchy' fluid mixing the commonly used Gassmann-Wood model (Gassmann, 1951) systematically underestimates partially saturated seismic velocities and often empirical models (Brie et al., 1995) are used to match the seismic data.
In an ideal situation either wireline sonic logs or laboratory experiments would be used to determine a velocity-saturation relationship that could be applied with confidence to seismic data at the fieldscale. However at Sleipner there is a lack of both wireline logs (shear wave velocities, for example, are limited to the Norwegian well 15/9-A23) and experimental data from multiphase flooding tests on Utsira cores in the laboratory (Arts et al., 2004). Such information is crucial to build up the detailed reservoir models necessary to robustly quantify in situ stored CO 2 .
In this paper, we present the first coupled geophysical and geomechanical measurements on an Utsira Sand core sample during a CO 2brine flooding experiment. We measure longitudinal (primary) acoustic wave velocity and attenuation, and electrical resistivity together with a continuous record of axial and radial strains. The test covers seven drainage brine-CO 2 fractional flow episodes simulating the transitional stages of increasing CO 2 saturation in a storage reservoir, plus an additional forced imbibition (brine flow), after cessation of CO 2 injection. For each flooding episode, pore pressure variations are induced in the system to investigate the geomechanical effect of inflation/depletion scenarios. We also use electrical resistivity as an indicator of CO 2 saturation and residual trapping. Our geophysical results are extrapolated to seismic frequencies using the recent rock physics model of Papageorgiou and Chapman (2017). The model simultaneously fits the velocity and attenuation data accounting for both patchy saturation effects and squirt flow dispersion, and provides a velocity-saturation relationship that can be applied to the seismic frequency band that lies between the patchy and uniform saturation bounds.

Experimental setup
The experiment was performed using the laboratory setup described in Falcon-Suarez et al. (2016). The rig (Fig. 1a) is capable of reproducing reservoir conditions up to 65 MPa of confining and pore pressure, and temperatures up to 50 C, on 5 cm diameter sample plugs. The setup is configured to simultaneously measure ultrasonic wave velocity, wave attenuation and electrical resistivity, together with axial and radial strains, during the co-injection of up to two pore fluids under controlled flow rates.
The sample is housed in a triaxial cell core holder. Inside the vessel a 6 mm thick hydrogenated nitrile butadiene rubber (HNBR) sleeve isolates the core plug from the confining fluid. The sleeve was perforated by 16 stainless steel electrodes for electrical resistivity measurements, radially distributed in two rings of 8 electrodes around the plug, and connected to an electrical resistivity tomography data acquisition system (ERT, Fig. 1a). Under typical operating conditions the resistivity measurement error is ∼5% (at frequencies 1-500 Hz) for samples in the electrical resistivity range 1-100 Ω m (North et al., 2013).
The sample is axially confined by two platen-PEEK buffer rod sets, which house ultrasonic pulse-echo sensors and implement fluid pathways that allow pore fluid circulation across the sample. These buffer rods have well defined acoustic impedance and low energy loss, providing a reliable delay path to enable the identification of top/base sample reflections for calculating ultrasonic P-and S-wave velocities (V p and V s ) and attenuations (inverse quality factors Q p −1 and Q s −1 ) using the pulse-echo technique (Best et al., 2007;McCann and Sothcott, 1992). The technique provides useable frequencies between 300 and 1000 kHz with absolute accuracies of ± 0.3% for velocity and ± 0.1 dB cm −1 for attenuation (Best, 1992). A hydrostatic confining pressure configuration was adopted for the experiment with both confining and pore pressure controlled by dual ISCO EX-100D systems; an extra ISCO ED100 cylinder being used as backpressure downstream of the sample. Pore fluid is indirectly delivered/received using fluid transfer vessels (FTVs, Fig. 1a) to avoid the contact between the controllers and brine/CO 2 fluids (i.e. minimizing corrosion effects), and additionally monitored with pore pressure (piezoresistivity) sensors located upstream and downstream of the sample. The FTVs are immersed in a thermal bath to set pore fluid temperature at the target conditions. For this experiment, the delivery FTVs are filled with pure CO 2 and brine (NaCl brine or CO 2 saturated NaCl brine depending on the experimental stage; see experimental procedure below), whereas the receiver vessel contains the resulting fluid. Two additional external pressure vessels store the injective fluids for refilling the FTVs (i.e., brine tank and CO 2 -brine vessels, Fig. 1a).

Sample preparation
The sample of Utsira Sand used in this study was cored from the upper section (E641 − 1085-1085.25 m) of the well 15/9-A-23, Viking Graben, Central North Sea (Table 1). The core sample was sandwiched by two annuli (PEEK) and two filter papers with pore size < 250 μm (Fig. 1b). This special preparation was adopted to counteract the very low consolidation state of the sample in order to (i) minimize mechanical damage during the assembly of the experiment, (ii) ensure an appropriate repartition of the axial loading and (iii) mitigate potential sand grain migration and pipe clogging effects. A 90°biaxial strain gauges set was epoxy-glued from one annulus to the other, for measuring axial and radial strains during the experiment. Sample mineralogy was determined by X-ray diffraction (XRD) analysis (Table 1).

Brine-CO 2 flow-through test
In order to simulate the progressive increase of the CO 2 to brine ratio in the storage reservoir during CO 2 injection, the experiment was set up as a two-phase flow steady state drainage test in which an increasing CO 2 to brine ratio flow is forced to pass through the sample (Falcon-Suarez et al., 2017;Falcon-Suarez et al., 2016;Perrin and Benson, 2010). For each CO 2 -brine co-injection episode, inflation/depletion cycles were simulated by stepwise variation of pore pressure. The test methodology and conditions replicate those used in Falcon-Suarez et al. (2016) and Falcon-Suarez et al. (2017) for simulating CO 2 storage, respectively in tight and weakly-cemented shallow siliciclastic reservoirs.
To replicate the state of stress at the top of the Utsira Sand at Sleipner, the lithostatic (overburden) stress was calculated at 16.4 MPa based on wireline overburden rock density measurements. In line with this, the experiment was performed at constant 16.4 MPa confining stress (σ 1 = σ 2 = σ 3 = σ c ). A scattering of fluid pressure measurements across the Utsira Sand broadly suggests initial hydrostatic conditions (Baklid et al., 1996). Pore fluid pressure (P p ) was initially set at 7 MPa, giving an initial effective pressure (P eff = σ c − P p ) of 9.4 MPa. The pore pressure was then increased in 1 MPa steps to 12 MPa, reducing effective pressure to 4.4 MPa, and then returned back to 7 MPa, completing the inflation/depletion stress-path. The geophysical parameters were measured on each step once the steady state condition was reached, based on the stabilization of the outlet flow record in the back pressure controller. The stress-path sequence was repeated eight times over seven consecutive drainage episodes followed by an imbibition one. During the drainage part of the test, brine:CO 2 fractional flow (Q brine :Q CO2 ) was progressively increased by 20% increments from 0 to 100% CO 2 , but keeping the total flow (i.e., Q brine + Q CO2 ) constant at 0.5 mL min −1 . The concluding imbibition stage comprised a brine flowthrough cycle after the last drainage episode.
Most of the Utsira CO 2 plume is at reservoir conditions above the Critical Point for CO 2 (31.1°C, 7.39 MPa), and for the experiment temperatures were kept above this value, oscillating between 31.1-32°C. Pressures were similarly maintained above the Critical Point except for the very first pore pressure setting which was intentionally lowered to 7 MPa to obtain measurements just below. Further details of the experimental procedure can be found in Falcon-Suarez et al. (2017).
During the test, axial and radial strains were measured continuously, every second. Electrical resistivity and ultrasonic measurements were systematically acquired at the end of each stress-path step, to obtain comparable values between both parameters.

Electrical resistivity into degree of saturation
Electrical resistivity was inverted using software based on the EIDORS MATLAB toolkit (Andy and William, 2006) for uniform/ homogeneous isotropic resistivity and heterogeneous isotropic resistivity distributions. Both inversion schemes employ a finite-element forward model of the sample and electrodes.
For our CO 2 -brine system, the degree of brine saturation (S w ) of the sample can be computed through the bulk electrical resistivity (R b ) using Archie's law (Archie, 1942), as follows: where R w is the electrical resistivity of the pore water, ϕ the porosity of the sample, and m, n and a are fitting parameters corresponding to the cementation factor, the saturation exponent and proportionality constant, respectively. Assuming that the system remains chemically invariable, and considering negligible the effect of the dissolved CO 2 on the electrical resistivity for such high salinity solution (Börner et al., 2013), then the above expression can be normalized with respect to the original (brine saturated flooding episode) bulk electrical resistivity of the sample (R 0 , when S w = 1) as follow (Carrigan et al., 2013): Pore fluid samples were collected downstream at different times during the experiment to study the effect of pore fluid changes on the electrical resistivity (Falcon-Suarez et al., 2017). The collected fluids were measured by inductively coupled plasma optical emission spectrometry.

Experimental results
The experiment lasted ∼8 days during which the flow-through test  was actively running ∼76 h (injection effective time), resulting in ∼120 pore volume (PV) throughputs. During this period the sample was subjected to 8 brine-CO 2 co-injection stages and 75 states of stress. The raw data obtained during the test (Fig. 2), show the seven drainage steps and the subsequent forced imbibition stage (R-100:0). For each brine:CO 2 fractional flow and stress-path step, V p and Q p −1 (obtained at 600 kHz from Fourier analysis of broadband signals) are plotted, together with the measured strains and associated porosity changes, and electrical resistivity. During the first two flow-through episodes (100:0 and 100 (s) :0, for brine and CO 2 -saturated brine, respectively) neither V p nor resistivity showed detectable changes due to the presence of dissolved CO 2 (blue band in Fig. 2). V p decreased with effective pressure from 2090 m s −1 to around 2000 m s −1 , and seemed to recover fully as pore pressure was reduced; Q p −1 was less affected by effective pressure changes. The arrival of free-phase (i.e., non-dissolved) CO 2 (corroborated by the resistivity increase; transition from blue to yellow band) triggered a sharp drop in V p (by ∼12%) followed, with increasing CO 2 fractional flow, by a further gradual decrease to ∼27% below the initial velocity. This minimum is similar to results obtained by Lei and Xue (2009) for Tako Sandstone, and slightly lower than those reported by Alemu et al. The volumetric deformation, equivalent to porosity variation (Δε v = Δϕ), was calculated from the axial (ε a ) and radial (ε r ) strains using the relationship ε v = ε a + 2ε r . The deformation ε v reached a maximum of ∼0.4% at for CO 2 fractional flows of 0.4 and above, which agrees with the results of Lei and Xue (2009). Quasi-total strain recovery during the stress sequences suggests the sample behaved approximately elastically for all brine:CO 2 fractional flow episodes during the drainage part of the test, except for the first injection of CO 2 (80:20) when the sample experienced an inflation of ∼0.15%; partially recovered during final imbibition cycle. This observation will be discussed in the geomechanical section below.
Likewise, the electrical resistivity also increased with CO 2 , from ∼1.5 Ω m for brine, up to ∼3.8 Ω m with pure CO 2 flowing through the sample (brine:CO 2 0:100). The sharpest change, up by ∼30%, took place during the initial 80:20 brine:CO 2 drainage cycle. No differences were measured between pure brine and the CO 2 saturated brine flooding episodes. Additionally, a small but consistently jump upwards occurs every time the pore pressure is lowered to 7 MPa from above (at the end of a pressure cycle). It could be related to local brine displacement effects due to decompression of CO 2 during supercritical-togas phase change.
During forced imbibition (R-100:0; green band) the brine flooded back into the sample and partially refilled the pore space. As a result, V p increased whereas Q p −1 , strains and resistivity decreased towards the initial values prior to injecting CO 2 . The imbibition shows a transitional evolution during the pore pressure increase, when brine replaced most of the CO 2 , preferentially affecting V p, Q p −1 and resistivity. Brine-CO 2 flow-through test on sandstone sample from Utsira. Pore pressure (P p ), effective pressure (P eff ), temperature (T), ultrasonic P-wave velocity (V p ) and attenuation factor (Q p −1 ), axial and radial strains (ε a and ε r ) and volumetric strains (ε v ), porosity (ϕ) variation, and electrical resistivity for eight consecutive brine:CO 2 fractional flows, covering seven drainage (the first, 100:0, using brine as pore fluid and the next six, from 100 (s) :0 to 0:100, using CO 2 saturated brine) and a forced imbibition (R-100:0) episodes, plotted versus pore volume (PV). Dark striped bands are the interludes between two consecutive brine:CO 2 episodes. Blue and yellow bands indicate drainage measurements, prior to and in the presence of free (non-dissolved) CO 2 , respectively (dark blue for brine; light blue for CO 2 saturated brine), and green for imbibition. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) I. Falcon-Suarez et al. International Journal of Greenhouse Gas Control 68 (2018) 236-246 Subsequently the pore pressure drop led to significant decrease of strains and minor V p and resistivity variations.

Resistivity into degree of saturation
Pore fluid analysis showed very little variation in the electrical resistivity of pore water with major cations concentrations consistent within the range 558 ± 4.4 mmol L −1 (see detailed chemical composition in the Appendix A), so we assume that any measurement error will lie within the instrumental error for resistivity (5%). Furthermore, the actual compositional analysis of the pore water may have been altered by the contamination with drilling fluids during coring and, therefore, we are not undertaking interpretations regarding geochemical changes on the original rock sample.
The resistivity measurements were converted into brine/CO 2 saturation (i.e., S w = 1 -S CO2 ), using the simplified form of Archie's law (Eq. (2)). The observed resistivity range, from ∼1.5 Ω m to ∼3.8 Ω m, corresponds to a S w range of 1-0.6 ( Fig. 3) with a near-linear best-fit curve (n = 1.8). The derived S w values were then plotted against V p and Q p −1 (Fig. 4). Albeit with significant scatter, there is a clear positive relationship between V p and S w (Fig. 4a), roughly consistent with previous rock physics estimates of the V p -S CO2 relationship for Utsira Sand which indicated minimum V p values within the S CO2 range 20-100% (Arts et al., 2004). Q −1 also varies with saturation, being lowest with a single pore fluid (S w = 1). We estimate the maximum saturation of CO 2 (S CO2,m ) attained in the sample at the end of the drainage stage to be ∼38%. Conversely at the end of the final imbibition stage, the residual CO 2 saturation (S CO2,r ), relating to the maximum brine saturation post-imbibition, was estimated by wet-dry weight difference to be < 4%. On the other hand, the final resistivity measurements give an estimate of residual CO 2 saturation of ∼6%. The weight-derived value might be reduced by exsolution during pressure release while removing the sample, although Zuo et al. (2012) found that the mobility of disconnected CO 2 bubbles within the brine is very low for rapid pressure drops higher than ours (7 MPa).

Simultaneous V p and Q p −1 and fluid distribution modelling
In this section, we aim to model the variation of V p with saturation Fig. 3. P-wave velocity (V p ) and brine saturation (S w ) versus resistivity. To facilitate visualization, one unique cross error per data source is displayed. Fig. 4. The best fit curves for velocity and attenuation based on the full saturation data. The best fit curves for simultaneous inversion of ultrasonic velocity and attenuation data.
The pressure-dependence of the data is not taken into account and inversion is performed for the patch parameter assuming gaseous (7 MPa) and supercritical (12 MPa) CO 2 but a common crack density. The black curves are the fits to the ultrasonic data, the red curves are their low-frequency limits, and the Gassmann-Wood limits (q = 1) for each CO 2 state define the shaded yellow area. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)  (2015); K m and μ s , estimated from XRD. b Brine properties from Batzle and Wang (1992); CO 2 properties from Span and Wagner (1996). of CO 2 , which requires an adequate modelling description of both the measured (ultrasonic) frequencies and the seismic frequency limit. This will allow us to extrapolate the variation of ultrasonic V p with saturation (Fig. 2) into the seismic frequency band.
To this end, we use the squirt flow model introduced by Papageorgiou and Chapman (2017) that combines the effects of wave induced fluid flow, which is dominant over the ultrasonic wave frequencies of the experiment (Amalokwu et al., 2015;Amalokwu et al., 2017), with patchy saturation (e.g., White, 1975), which may affect the seismic wave propagation through the sample.
One of the simplifying assumptions we make is that the ∼5% variation in the fully brine saturated velocity observed between extremes of the pressure cycle is averaged over. This means that we will not account for pressure-induced changes in the frame moduli in this section.
Another assumption concerns the dry moduli of the Utsira sample. As it was not possible to measure velocities in the dry sample, we take values from the literature ( Table 2) chosen in accordance with the velocities observed in the Utsira reservoir by Furre et al. (2015). Our aim is then to compare the seismic frequency limit of our model with their findings and provide a rock-physics based method for seismically quantifying CO 2 .
Even though the impact of the pore pressure on the matrix properties is averaged, the varying physical properties of CO 2 during the pressure cycles are significant and cannot be ignored. We therefore use the pore pressure endpoints at 7 MPa and 12 MPa as indicators of gas and supercritical CO 2 respectively and invert the model independently for each of these two cases. The physical properties of CO 2 at these pressures (and for the given temperature of the experiment) are obtained from the NIST database and are shown in Table 2 together with the brine properties under the same conditions. With these inputs for the fluids, we use the model of Papageorgiou and Chapman (2017) to calculate an effective fluid modulus and characteristic squirt flow frequency. These values depend on the bulk moduli and viscosities of the two saturating fluids, which means that different states of CO 2 produce different effective fluids and frequencydependence. Explicitly, the effective fluid modulus K f and characteristic frequency ω c of this theory are given by: (3) and = ⎛ ⎝ and where K w and K CO2 are the respective brine and CO 2 moduli, k w and k CO2 the relative permeability values (from Papageorgiou and Chapman, 2017) for the flow of brine and CO 2 through the Utsira sand, and η w and η CO2 the respective viscosities of the two fluids. Here we have scaled ω c against the characteristic frequency ω 0 of brine, assumed equal to the experiment frequency (600KHz). This choice is driven by observational evidence of past work in sandstones (Amalokwu et al., 2017), where the frequency band for the transition of the modulus from its relaxed to unrelaxed state is observed around the experiment frequency, but arbitrary in its essence. More important are the relative changes in the characteristic frequency when the matrix is partially saturated. Explicitly, depending on the patch parameter q, the unrelaxed state may be observed under partial saturation whereas the fully saturated matrix is still in a transition.
The parameter q is a bounded constant free parameter representing the averaged departure from pressure equilibration across the pore space and can be justified by considering patch-style effects, capillary effects, membrane stresses, etc. Mathematically, it is expressed via: where ΔP w , ΔP CO2 represent the averaged fluid pressure response to stress. In the bound q = 1 the fluids respond uniformly to stress and the mixing law of Eq. (3) corresponds to Wood's law, whereas for smaller values of q, the effective modulus in equation 3 resembles Brie's law (see Papageorgiou et al., 2016).
With the above definitions, one can calculate P-wave seismic velocity (V p ) and attenuation (Q p −1 ) using: where ρ is bulk density, μ d is the shear frame modulus and i is the imaginary unit. The frequency-dependent K eff is read from Papageorgiou and Chapman (2017): and In these equations, φ p , φ c are the pore and crack porosities, ε and r are the crack density and aspect ratio, respectively, and K m , μ and ν are the grain bulk and shear moduli and Poisson's ratio respectively. The first two lines in the expression for the bulk modulus are independent of ω and they sum to Gassmann's model with the effective fluid of Eq. (3). With these definitions we assume the crack density ε and patch parameter q to be fitting parameters to be inverted for using the ultrasonic data. The inversion is based on a non-linear Nelder-Mead optimization scheme. The velocity and attenuation data are inverted simultaneously. Also, because of unclear evidence of pore pressure dependency for the two parameters during partial CO 2 saturation stages (yellow band, Fig. 2), we performed the inversion twice over the entire dataset, disregarding pressure dependence in the data: once assuming the most compliant CO 2 (at 7 MPa) and once the stiffest CO 2 (at 12 MPa). In performing the inversion the crack density is assumed to be independent of fluid but we allow the parameter q to differ between inversions (Table 2). For each assumed CO 2 pressure, we obtain three different curves associated with the inversion (shown in Fig. 4): the ultrasonic velocity and attenuation curves fitting the data for 7 MPa and 12 MPa, their respective low frequency limits and the Gassmann-Wood (uniform saturation; upper and lower bounds enclosing the yellow shaded area in Fig. 4).
Our prediction for the seismic velocity variation with CO 2 saturation (red curves in Fig. 4), strongly agrees with the upper bound of Gassmann-Wood, which is reflected in the inversion result q = 0.79 for the inversion at 12 MPa (where q = 1 corresponds to the Gassmann-Wood limit). However, the lower bound of our prediction differs significantly from the lower bound of the yellow area (corresponding to q = 1, the Gassmann-Wood limit for CO 2 at 7 MPa), which indicates more patchy distribution for gasseous CO 2 (q = 0.09).
The common crack density (ε, Table 2) indicates that the amount of dispersion is relatively similar for either the gas or supercritical CO 2 state, which has implications for the quantification of CO 2 . For instance, Furre et al. (2015) noted that the observed time-shifts from Sleipner correspond to layer thickness < 30m. In turn, it implies that CO 2 saturated sand velocities can vary between 1.4 and 1.5 km s −1 . Based on our seismic velocity-saturation model, this observation constrains the CO 2 saturation in the Utsira sand to be between about 10% and 50% (where these two velocity bounds meet our red dashed and solid curves in Fig. 4). Other studies (Chadwick et al. (2016)) suggest lower velocity estimates which, based on our modelling study suggests CO 2 saturation > 20%.

Geomechanical assessment
The data collected during this experiment allow us to investigate the interplay between mechanical and pore fluid distribution effects, using the electrical resistivity as an indicator of pore fluid distribution (Fig. 5). For values of S w above 0.9, V p is a good geomechanical indicator, showing a clear linear relationship with strain (R 2 = 0.738); although the available data within this saturation range lie mostly below the estimated residual CO 2 saturation value of 6% (i.e., S w > 0.94). For S w within the range 0.8-0.9, small changes in the pore fluid composition lead to large V p variations (Fig. 4a), which mask mechanical effects; while, for S w < 0.7, sample deformation correlates only poorly with V p (R 2 = 0.241). On the other hand, Q p −1 shows significant data scattering for the whole saturation range.
To further investigate pressure effects in the Utsira Sand we focus on V p -ε v relationships. With the exception of the initial input of CO 2 (80:20 brine:CO 2 fractional flow episode) when there was a prompt geomechanical inflation (Fig. 2), and the reverse during the final forced imbibition, the sample behaved in an approximately elastic fashion, with no permanent deformation. The data therefore appear suitable for inferring pore pressure changes from V p and the associated ε v . The geomechanical inflation associated with the arrival of the CO 2 has been previously reported in CO 2 flooding experiments (Falcon-Suarez et al., 2017;Falcon-Suarez et al., 2016) and also, at field scale, at the In Salah storage site in Algeria (Mathieson et al., 2011;Onuma and Ohkawa, 2009;Vasco et al., 2008), where the observed geomechanical deformation is consistent with CO 2 injection induced-pressure built-up forward and inverse modelling (Mathieson et al., 2011). In our experiment, pore pressure increase alone is insufficient to explain the permanent deformation because it only occurs after CO 2 is present. The strain gauges are embedded in a carrier chemically inert to CO 2 (with no evidences of equipment damage after the test), whereas the conditions (pressure and temperature ranges) are repeated from brine to the CO 2 rich flooding cycles. Thus, in the absence of additional information, we assume the effect is not an experimental artefact. Instead, we interpret this phenomenon as the interplay between the pore pressure increase and CO 2 -induced salt (NaCl) precipitation, although further investigation is recommended in this regard.
Complex salt precipitation patterns associated with the supersaturation of the original brine due to evaporation of water into the CO 2 stream have been reported by a number of studies, including experimental works (Miri et al., 2015) and numerical simulations (Kim et al., 2012). This phenomenon has been found to be a rapid process that preferentially occurs around the main flow channels (Ott et al., 2011), Fig. 5. (a) P-wave velocity (V p ) and (b) P-wave attenuation (Q p −1 ) versus volumetric deformation (ε v ) every 10% brine saturation (S w ) intervals. Linear fittings were obtained for S w > 0.9 (L Sw > 0.9 ) and S w < 0.7 (L Sw < 0.7).

Table 3
Fitting parameters for the expression proposed by Eberhart-Phillips et al. (1989) to relate P-wave velocities with effective pressure (Eq. (14)). leading to significant porosity and permeability reduction (Bacci et al., 2013;Jeddizahed and Rostami, 2016). Furthermore, salt crystallization pressure can locally reach values even over 150 MPa in confinement (Desarnaud et al., 2016). When the salt crystallization pressure exceeds the pore pressure, the effective pressure condition of the rock is reduced (Zheng et al., 2015), resulting in some volumetric dilation. In our sample, this process would be enhanced by the pore pressure increase path (from 7 to 12 MPa) during the 80:20 brine:CO 2 fractional flow episode. Induced pore dilation would allow salt precipitation between sand grains, resulting in permanent deformation. Afterwards, the dilated sample behaves elastically until the imbibition episode, when brine dissolves salt aggregates and the sample approximately recovers its original volume. Prior to this work, no geomechanical experimental data were available for the Utsira Sand at Sleipner, so Chadwick et al. (2012) applied the empirical expressions proposed by Eberhart-Phillips et al. (1989), relating seismic velocities (V p ) to effective pressure, i.e., P eff for Utsira sand, within the fully water −saturated part of the reservoir (S w = 1). Such expressions were obtained from an experimental dataset of water-saturated sandstones, and depend on the porosity (ϕ) and clay content (C). The original form of the equation is: p e f f gP eff (14) with a = 5.77, b = 6.94, d = 1.73, f = 0.446 and g = 16.7 the fitting parameters for compressional waves (Eberhart-Phillips et al., 1989). Using the clay content obtained by XRD (C = 0.07; Table 1), and the porosity estimated for our sample (ϕ = 0.375), we adapted the above expression to our experimental data using nonlinear regression, firstly for the brine-saturated sample and then for a range of brine saturations at 10% intervals (Table 3, Fig. 6a). For the brine saturation condition, measured and calculated parameters strongly correlate (R 2 > 0.85), giving a V p :P p rate of ∼14.6 ± 0.6 m s −1 MPa −1 for the considered P eff range. This value is approximately 50% lower than that derived by Chadwick et al. (2012) from Eq. (14).
In addition, to quantify the geomechanical effect we have adjusted the deformation suffered by the sample due to pore pressure, using linear regression (Fig. 6b). Note that we assume effective pressure coefficient equals to one, the accepted value for very weakly-cemented sandstones (Hofmann et al., 2005).

Discussion
The measured value of V p for the brine-saturated samples was around 2090 m s −1 at 7 MPa (Fig. 2). This is compatible with in situ V p measurements from wireline logs in the Utsira Sand, which range from ∼1950 to ∼2100 m s −1 with a quoted average of 2056 m s −1 (Zweigel et al., 2000). This fact suggests that the core sample is reasonably representative of the Utsira Sand as a whole, but also that either no significant frequency effects take place, or several dispersion mechanisms cancel each other.

Pore fluid distribution
We have analysed the effects of both fluid distributions and pore pressure changes from the experimental data. To account for patchy pore fluid distribution, we simultaneously inverted the ultrasonic Pwave velocity and attenuation data using a rock physics model based on combined squirt flow and patchy saturation effects − both known to be significant factors affecting the dispersion of partially CO 2 saturated sandstones from previous experimental studies (e.g., Alemu et al., 2013;Falcon-Suarez et al., 2016). The low frequency V p −saturation relationship derived from the ultrasonic inversion depends on the physical state of CO 2 (whether it is below or above the Critical Point). In either case, however, the seismic velocity-saturation relation derived from the inversion is bounded by the two red curves in Fig. 4. It should be noted that we also performed the inversion assuming intermediate values of CO 2 pressure and the inverted curves always fell within those bounds so we chose to present only the extreme cases corresponding to gasseous (7 MPa) and supercritical (12 MPa) CO 2 . The inversion for the parameter q seems to largely compensate the compliance of the soft CO 2 , where for 7 MPa it inverts to q = 0.09, which corresponds to a relatively patchy signature. However, this effect is exaggerated by ultrasonic dispersion (black dashed cuve in Fig. 4). As a result, for seismic frequencies where no squirt flow dispersion is expected, our predicted lower limit for the velocity-saturation relation (red dashed cuve in Fig. 4), falls between the Gassmann-Wood prediction for the two states.
In this regard, several works have aimed at quantifying in situ CO 2 stored in the Utsira Sand using a variety of methods based on time-lapse 3D seismics (e.g., Chadwick et al., 2005;Furre et al., 2015;Williams and Chadwick, 2012). The CO 2 plume at Sleipner has a tiered structure Fig. 6. (a) P-wave velocity (V p ) and (b) volumetric deformation (ε v ) versus pore pressure increment (ΔP p ) every 10% brine saturation (S w ) intervals and for the brine saturated condition (S w = 1).
I. Falcon-Suarez et al. International Journal of Greenhouse Gas Control 68 (2018) 236-246 comprising a number of sub-horizontal reflections interpreted as thin layers of CO 2 trapped at various levels within the reservoir (Arts et al., 2004). A key uncertainty in quantification is the state of CO 2 saturation within these layers. This can be surmised from laboratory-determined capillary pressure measurements on core samples (Chadwick et al., 2005), but plume-scale velocity determination directly from the 3D seismics would likely be more representative and diagnostic of plume flow processes. This however is very challenging. A number of seismic inversion approaches have been tested on the data (e.g., Clochard et al. (2009)), but inability to resolve the thin layers properly has always compromised velocity determinations. Chadwick et al. (2016) used an interpretive inversion scheme based on structural analysis of the CO 2 -water contact to estimate V p for the topmost layer in the range 1350-1430 m s −1 , with an uncertainty of the order of ± 100m s −1 . At hydrostatic conditions , pressure at the reservoir top would be around 8 PMa, so reference to the low frequency curves (Fig. 4a) indicates that the velocity is consistent with CO 2 saturations in the range 0.2-1. This rather wide range is consistent with previous laboratory capillary pressure tests on Utsira Sand core which indicate residual water saturations as low as 0.05 (Chadwick et al., 2005), but also with our own test which suggests residual water saturations of 0.06. The reason for this estimated residual water saturation discrepancy is unclear, but might represent real differences in core samples, or perhaps experimental method.
A CSEM survey has also been deployed In Sleipner, but technical and data-processing issues have complicated the interpretation (Park et al., 2013;Park et al., 2017). The resistivity data collected during our experiment might improve qualitative and quantitative determination of the CO 2 distribution in the aquifer from electromagnetic modelling (Park et al., 2017).

Pore pressure and geomechanical effects
The Utsira Sand is interpreted as a large hydraulically well-connected aquifer with pressure increase at the injection well up to 2006 reported as < 0.2 MPa  with very little geomechanical deformation at the current injection rate (Verdon et al., 2013). However, thin low permeability barriers (intra-reservoir mudstones) have been recognized from wireline logs in the vicinity and also from the time-lapse seismics, so local pressure increases related to unexpected compartmentalization cannot be ruled out. Chadwick et al. (2012) used statistical analysis of very small time-lapse time-shifts to constrain pore pressure variations from V p , and following this work we moved forward to calibrate the empirical relationship proposed by Eberhart-Phillips et al. (1989) but taking into account fluid saturation as well (see above). The best fits between V p and pore pressure were found at either very low or high CO 2 saturations, while values in between were masked by the stronger effect of pore fluid saturations.
We also analysed the mechanical deformation related to pore pressure variations in the sample. We found that the linear slopes Δε v / ΔP p are practically constant with a slight tendency to decrease with increasing CO 2 saturation, but this effect might be related to sand grain reorganisations in the sample after the initial inflation/depletion cycles, as suggested by the reduction in sample deformation from the first to the second brine flow episode (Fig. 2). By contrary, CO 2 diffusion within clay (particularly kaolinite) layers may contrarily affect the rock stiffness of the rock, by altering (weakening) the microstructure of clay particles (Delle Piane and Sarout, 2016). However, this effect can be also explained by the CO 2 -induced salt precipitation hypothesis, which supports the strain dependencies observed in Figs. 5 and 6, slightly stronger for high (L Sw > 0.9) than for low (L Sw < 0.7) brine saturations. In essence, advanced CO 2 saturations states led to a salt-skeleton reinforced rock frame, increasing the rock stiffening.
The slope Δε v /ΔP p also provides information about the 'pore compressibility' of the porous medium, defined as B p = ϕ −1 (Δε v /ΔP p ). B p is normally a poorly constrained parameter for most of reservoirs formations due to its high variability Cui et al., 2010), but of great importance in coupled geomechanical and flow modelling to predict the elastic deformation of weakly consolidated rocks, fluid pressure increase and related porosity-permeability variations (Minkoff et al., 2003). Our experimental data give a value of B p ∼1.2 × 10 −10 Pa −1 for the Utsira Sand, based on the relationship for brine saturation. This value is similar to those published for this parameter in generic reservoir rocks (Freeze and Cherry, 1979); although our estimate might be lower than the true value because (i) the experimental determination relies on the assumption of hydrostatic confining stress and (ii) our laboratory setup combines small sample and high sensibility pumping controllers, which rapidly dissipates the confining (pore pressure-induced) overstress. Verdon et al. (2013) point out the importance of a comprehensive geomechanical assessment of the target reservoir before CO 2 injection, remarking the importance of laboratory measurements. The significance of our results is in combining real geophysical and geomechanical properties of the Sleipner storage site, including baseline data (prior to CO 2 injection).

Dissolution effects
The brine and CO 2 -saturated brine flow-through tests indicate that dissolved CO 2 , up to 9 vol% at the experimental conditions (Berg et al., 2013), is not detectable seismically or by resistivity, confirming previous published assumptions in this regard (Börner et al., 2013;Chadwick et al., 2010;Eiken et al., 2011). The lack of a resistivity signature is compatible with the known salinity in the Utsira Sand and contrasts with the situation at the Nagaoka pilot site in Japan, where time-lapse resistivity logging has been able to track the gradual dissolution of an injected CO 2 plume, because the reservoir has very low salinity (Mito and Xue, 2011).

Conclusions
An experimental geophysical and geomechanical investigation of the Sleipner CO 2 storage reservoir has been carried out for the first time. Working with the very weakly-cemented Utsira Sand at reservoir conditions is experimentally very challenging and previous attempts have not been successful. The resulting dataset combines measurements from variable brine-CO 2 saturation and pore pressure conditions. From these we have (i) calibrated a rock physics model and developed a velocity-CO 2 saturation relationship that can be applied to the seismic band that lies between the patchy and uniform mixing bounds, and (ii) inferred empirical relations between pore pressure changes and geomechanical effects in the reservoir, for different saturation ranges.
Our experimental results represent therefore a unique source of data to complement the existing geophysical field data (seismic and electrical resistivity) and geomechanical information on Sleipner CO 2 storage site.