Impact of Fluid Distribution and Petrophysics on Geophysical Signatures of CO2 Storage Sandstone Reservoirs

Carbon Capture and Storage (CCS) is a key element to achieving net‐zero energy challenge timely. CCS operations require the integration of geophysical data, such as seismic and electromagnetic surveys, numerical reservoir models and fluid flow simulations. However, the 10–100s m resolution of seismic imaging methods complicates the mapping of smaller scale rock heterogeneities, while borehole measurements commonly show large fluctuations at sub‐cm scales. In this study, we combine laboratory data, well‐logging, rock physics theories and a proof‐of‐concept time‐lapse seismic modeling to assess the effect of pore‐scale fluid distribution and petrophysical heterogeneities on the expected performance of whole‐reservoir CCS operations in deep saline aquifers, by analogy to the Aurora CCS site, North Sea. We monitored the elastic and electrical properties of three sandstone samples with slightly different physical and petrographic properties during carbon dioxide (CO2) flow‐through tests under equivalent in situ effective pressure. We inferred the CO2‐induced damage in the rocks from the variations of their hydromechanical properties. We found that the clay fraction, CO2‐clay chemical interactions, and porosity were the main factors affecting both the CO2 distribution in the samples and the hydromechanical response. We used seismic modeling of well‐log data and the laboratory results to estimate the reservoir‐scale time‐lapse seismic response to CO2 injection and to assess the effect of the rock heterogeneities in our interpretation. The results show that disregarding the effect of rock heterogeneities on the CO2‐brine fluids distribution can lead to significant misinterpretations of seismic monitoring surveys during CCS operations in terms of both CO2 quantification and distribution.


Introduction
Net zero carbon emission technologies are essential to meet Paris agreement targets (UNFCCC, 2015), including offshore energy storage (associated with the hydrogen fuel economy) and Carbon Capture and Storage (CCS) schemes; the former is still a promising technology that will benefit from the established expertise on CCS.In Europe, the oil and gas industry majors commitment to develop CCS infrastructure in both ongoing and planned offshore carbon dioxide (CO 2 ) storage in the North Sea is crucial to enable the energy transition (Global-CCS-Institute, 2021).
Norway is a pioneer country in offshore CO 2 storage.In 1996, Sleipner Field, central North Sea, became the first commercial CCS project in history, with an injecting rate of ∼1 MTn CO 2 per year.Up until now, more than 20 MTn of CO 2 has already been injected into the Utsira sandstone reservoir (Dean et al., 2020).Sleipner was followed by the Snohvit project in 2008, with an injection rate of ∼0.7 MTn CO 2 per year.More recently, the Norwegian government, in partnership with oil and gas industry, has committed to large-scale CO 2 sequestration in the Norwegian Continental Shelf (Global-CCS-Institute, 2021).Promising sites have been identified in Smeaheia (Fawad et al., 2021) and Aurora (Aker et al., 2021) areas, with high-quality sandstone reservoirs with a certain degree of heterogeneity (i.e., variable clay content and porosity).
Assessing the impact of reservoir heterogeneities is a key factor in CO 2 storage reservoir predictions, as vertical and horizontal heterogeneities condition the injectivity and CO 2 plume migration paths (Sundal et al., 2013).Al-Khdheeawi et al. (2017) remarked on the impact of reservoir heterogeneity on the CO 2 migration and trapping capacities, while Williams and Chadwick (2021) have recently shown how reservoir-scale heterogeneities have conditioned the CO 2 plume distribution at the Sleipner Field.However, most reservoir heterogeneities are undetectable from standard geophysical remote sensing technologies used in reservoir exploration -seismic and electromagnetic surveys, in particular.The small size, continuity or shape of these heterogeneities complicate the interpretation of the injection activities from the geophysical record (Chadwick et al., 2004).
Reservoir heterogeneities are ultimately linked to the physical, mineralogical and textural properties of reservoir rocks, which also condition the CO 2 storage capacity and fluid(s) distribution (e.g., Muñoz-Ibáñez et al., 2019).Original clay content, grain size distribution, type and degree of cementation, porosity and permeability are among the most relevant parameters, particularly for siliciclastic formations that are commonly barely reactive to CO 2 .Clay content might rapidly vary between adjacent rocks, potentially affecting porosity and permeability (e.g., pore clogging), and our interpretation of both parameters from elastic and electrical data.For instance, rock density increases with clay content (and so P-and S-wave velocities), while excess of ions lowers the bulk resistivity even for sea-like salinity environments (Han et al., 2011).Also, the clay composition, content and distribution in the pores and rock frame affect the absorption of CO 2 (Hamza et al., 2021).
Laboratory testing using reservoir samples provides valuable data sets that complement information from well logs and field geophysical surveys.Laboratory tests use accurate tools to identify and interpret the geophysical signatures corresponding to CO 2 plume migration from controlled experiments under in situ reservoir conditions (e.g., Chadwick et al., 2019).However, core-scale laboratory experiments are barely representative of the events occurring at the field spatial and temporal scales.The short lifetime of CCS reservoirs with respect to other geological processes reduces the temporal scale, particularly in reservoirs with low chemical reactivity to CO 2 , such as saline siliciclastic aquifers (i.e., North Sea-like reservoirs).By contrast, spatial uncertainties remain, involving reservoir heterogeneity and anisotropy, lateral extension and connectivity between different geological units, and geophysical techniques used for monitoring purposes.
Here we study the impact of rock heterogeneities on CO 2 distribution and the CO 2 -induced hydromechanical processes that might occur in siliciclastic CCS reservoirs, at the laboratory scale.The paper is structured, as follows: First, we perform CO 2 flow-through tests in the laboratory using three sandstones with similar mineralogy but different physical, elastic and hydraulic properties.During these tests, we analyze the variations of their geophysical signatures (ultrasonic waves and electrical resistivity) with increasing CO 2 content under the specific effective stress conditions of the target CCS reservoir at Aurora (Aker et al., 2021).We discuss the challenges posed by temperature in rock physics laboratories and propose a simple method to accommodate real PT conditions of deep CCS reservoirs into experimental procedures at constant (room) temperature, using common rock physics methods and assumptions.We complement our experimental study with a combined elastic and hydromechanical assessment performed on the three samples before and after the CO 2 flow-through tests.
Second, we perform a proof-of-concept time-lapse seismic modeling to analyze how the presence of heterogeneities might condition our geophysical interpretation of CCS reservoirs.In essence, the three sandstones analyzed in this study are representative of different formations that make up a vertically heterogeneous reservoir.We conduct a 4D seismic feasibility test to assess the sensitivity of seismic methods to detect changes in these formations linked to CO 2 injection.We compare our results to the simpler case, that of a vertically homogeneous reservoir comprised of a single formation unit.This comparison highlights the importance of laboratory testing design to obtain high-quality data sets representative of reservoir conditions, and the further need to properly upscale the results using the correct rock physics assumptions.

CO 2 Storage Context: Aurora Site
The CO 2 storage site Aurora is a seawater brine-bearing sandstone reservoir, which comprises the Johansen and Cook Formations from the early Jurassic Dunlin Group.Both formations are interconnected through faults (Aker et al., 2021), covering a ∼120 m thick effective storage reservoir.The reservoir formations contain laterally continuous calcite cemented layers of low permeability (Sundal et al., 2013).The injection site (well 31/5-7) is located southwest of the Troll field and is currently being developed by Equinor in partnership with Total and Shell as part of the Northern Light project (Aker et al., 2021).Figure 1 shows the distribution of the properties of the reservoir formations from well-log geophysical data (well 31/5-7), and the likely storage formations of Aurora site framed in yellow (on which we will focus our study).
The Johansen Formation is a medium to fine-grained, micaceous, moderately-sorted (upper part) to well-sorted sandstone, with thin calcite cemented streaks.The Cook Formation comprises very fine to fine-grained, subangular to sub-rounded, well-sorted, hard to friable sandstones with variable content of mica, glauconite and carbonates, and calcareous cement (Vollset & Doré, 1984).Estimates from core analysis revealed that porosity ranges within 16%-25% and vertical permeability is below 25 mD, with a vertical-to-horizontal permeability ratio <0.2 near the well injection point (below 2500 mbsf) depth (Sundal et al., 2013).The geophysical characterization of the area was performed by combining wellbore and 3D seismic data (Aker et al., 2021;Fawad & Mondol, 2018).The resulting rock physics analysis revealed that (a) the lower limit of the porosity estimates from core analysis decreases down to 5%, (b) the clay-sized particles volume fraction (C) ranges from 0.1 to ∼0.4,(c) the reservoir contains calcite laminations and cementation, and (d) quartz cementation increases with depth (Fawad & Mondol, 2018).In this study, we focus on the heterogeneity effects induced by variable clay-sized fraction within the reservoir formation on the geophysical interpretation of the CO 2 distribution and the changes in its original elastic and transport properties.Well-log geophysical information from well 31/5-7: P-and S-wave sonic velocities (V P and V S ), bulk density (ρ), porosity (ϕ) and clay-sized particles volume fraction (C) derived from gamma ray (GR) normalized with respect to the considered GR range for each value (i) as C = (GR i GR min )/(GR max GR min ).Solid dots are average values of each parameter within a given geological formation.Yellow frame shows the reservoir section considered for modeling in this study (see main text).

Core Samples Information
We selected our samples from the rock repository at the Rock Physics laboratory at the National Oceanography Centre (NOC RPL), Southampton, considering the petrographic, physical, elastic and transport properties described for the Johansen and Cook Formations (see above), that is, the target formations of Aurora's CCS reservoir.We found rock samples with similar properties within the Forties Formation (sandstone), well 22/14a-2 in the Central Graben of the UK Continental Shelf, North Sea (Couch et al., 2020), around 2640 ± 30 mbsf.From there, we selected three samples (here referred to ST1, ST2, and ST3) with a volume of clay-sized particles (C, which includes both major and clay minerals) ranging from 10% to 26% and porosities above 0.2 (Table 1; Figure 2), to study the effect of both parameters on the elastic and transport properties during and after CO 2 injection.Sample mineralogy was determined by X-ray diffraction (XRD), porosity by He-pycnometry from the cores dried (65°C oven) before (ϕ 0 ), and after (ϕ 1 ) the CO 2 flow-through tests.
The clay-sized fraction volume C was estimated from gamma ray well-logging data from the UK National Data Repository https://ndr.nstauthority.co.uk/, and then corroborated by thin section digital image analysis using the free software Fiji/ImageJ.We applied the machine learning tool Trainable Weka Segmentation (Arganda-Carreras et al., 2017) to estimate grains, clay and pores, and obtained results in agreement with well log data and lab measurements (Figure 2).Also, we applied a median filter to isolate grains and analyze particle sizes.
The results of the thin section image analysis reveal that the distribution of the clay-sized particles and porosity are in good agreement with those estimated from well-logs and He-pycnometry (see Table 1).ST1 shows the largest deviation in porosity (∼5% higher than measured), which can be related to a more heterogeneous grainsize distribution at the mm-scale.

Brine and CO 2 Flow-Through Tests: Experimental Setup
The tests were conducted using the high-pressure room-temperature (20°C) experimental setup for multi-flowthrough tests in the NOC RPL, Southampton (Figure 3).The setup allows monitoring geophysical properties (i.e., ultrasonic waves and electrical resistivity) of 5 cm diameter rock (core) samples, subjected to the injection of different pore fluids.For each sample, we performed a CO 2 flow-through test with geophysical and hydromechanical monitoring, together with an assessment of their elastic and transport properties pre-and post-CO 2 injection.To satisfy the needs for the CO 2 flow-through experiments, the rig was set up in a multi-flow configuration under well-controlled flow confining and pore pressure conditions (ISCO pump controllers).For all tests, we adopted hydrostatic confining stress (i.e., σ 1 = σ 2 = σ 3 = P c ). Pore pressure (P p ) was supplied using fluid transfer vessels (FTVs) to minimize fluid-induced corrosiveness effects on the equipment.The rig implements sensors for measuring, simultaneously, ultrasonic P-and two orthogonally polarized S-waves attributes (velocity and attenuation), electrical resistivity and axial strains.The ultrasonic sensors are housed in the two platens that axially confine the sample through a polyether ether ketone (PEEK) buffer rods of welldefined acoustic impedance and low energy loss, which allow calculation of ultrasonic velocities (V P and V S ) to a precision of ±0.1% and accuracy of ±0.3% (95% confidence, respectively) and P-and S-wave attenuation coefficient to an accuracy of ±0.2 dB/cm giving attenuation expressed as the inverse quality factors (Q P 1 and ) of ∼5% at 600 kHz, the frequency used to generate the ultrasonic data in this study, obtained from Fourier analysis of broadband signals using the pulse-echo technique (Best, 1992).From the axial platens, two arms hold linear variable differential transformer (LVDT) sensors for axial strain monitoring.Inside the triaxial vessel, the rubber sleeve that isolates the rock sample from the confining mineral oil implements an array of 16 stainless steel electrodes.Once the confining pressure increases and the sleeve is in contact with the sample, we can measure the bulk electrical resistivity using an electrical resistivity tomography (ERT ) data acquisition system designed and developed at the NOC (North et al., 2013).Under our experimental P-T-fluid salinity conditions (i.e., ∼3.5% NaCl brine), bulk resistivity error is <1% for the (homogenous and isotropic) fully saturated porous sandstones, raising up to ∼5% when the fluid is heterogeneously distributed (Falcon-Suarez et al., 2021).We refer to North

Experimental Workflow
We performed an elastic and transport assessment on each sample (i.e., ST1, 2 and 3), before and after being subjected to CO 2 flow-through to investigate potential CO 2 -induced changes.This assessment included the measurement of ultrasonic V P and V S , and their attenuations Q P 1 and Q S 1 on the samples both dry and brine saturated, together with the permeability to N 2 when dry, and electrical resistivity when brine saturated.The analysis focused on the stress dependency trends (and the end-point absolute values) exhibited by every parameter through an effective stress (P eff = P c nP p , with n being effective stress coefficient n (Biot, 1962)) path 5-10-15-20-26 MPa under drainage conditions.We used a gentle loading/unloading stress rate of ∼0.01 MPa/s and waited for one hour for every P eff step to ensure the sample was in stress-strain equilibrium before taking measurements.
The CO 2 flow-through tests each combined four consecutive flow-through stages consisting of (a) brine flow, (b) CO 2 -saturated brine flow, (c) CO 2 flow-through the brine saturated sample (i.e., drainage stage simulating the CO 2 injection into the saline aquifer), and (d) brine flow back through the sample partially saturated in CO 2 (i.e., imbibition stage simulating the natural aquifer recharge).During the tests, the monitoring included ultrasonic wave attributes (V P , V S , Q P 1 , and Q s 1 ), electrical resistivity and axial strain.

Modeling Approaches
To explain our ultrasonic measurements, we use Gassmann's equation (Gassmann, 1951) with an effective fluid assumed at iso-stress conditions based on Wood's law (e.g., Falcon-Suarez et al., 2016), and the White and Dutta-Odé model (see Mavko et al., 2009) to model the effects of fluid patches.Thus, we are able to (a) estimate P-and S-wave velocities over the whole brine saturation range, at a frequency of 600 kHz, (b) investigate sub-core scale heterogeneities using appropriate averaging for the mineralogy, and (c) assess the temperature effect in the predictions, using accurate fluid properties from their equations of state.Fluid properties (i.e., density, viscosity and compressibility at the experimental pressure, temperature and salinity conditions) were calculated using the High-Pressure International EOS for brine (Millero et al., 1980), and by interpolating data from the National Institute of Standards and Technology (https://webbook.nist.gov/chemistry/fluid/)for CO 2 .For the mineral moduli, we used Voigt-Reuss-Hill averages (e.g., Mavko et al., 2009), based on the mineralogical and original porosity (i.e., prior to CO 2 injection) estimates displayed in Table 1.
The transformation of bulk electrical resistivity (R b ) into degree of brine saturation (S w ) was performed using the modified Archie's empirical relationship to account for the contribution of clay minerals, based on the Waxman-Smits-Juhasz model (Juhász, 1981).This model accounts for the excess of charge through the cation exchange capacity (CEC) and the increased mobility of ions along the clay surfaces when in contact with an electrolytic solution.We used the CEC values proposed by Thomas (1976) for the clay minerals presented in our samples (i.e., illite = 0.1 meq g 1 , kaolinite = 0.03 meq g 1 and chlorite = 0.01 meq g 1 ), and applied the Waxman-Smits-Juhasz model formulation as described in Mavko et al. (2009) and Falcon-Suarez et al. ( 2021) to obtain S w , as follows: with the resistivity formation factor F = aϕ m , m and n are the cementation and saturation exponents, a is an empirical parameter close to unity, and R w is the resistivity of the brine solution; for temperature (T ) in degrees Celsius, and where ρ s is mineral grain density.
Due to the theoretical low CO 2 -reactivity of the major minerals that form our samples (i.e., siliciclastic sandstones with no carbonates), we anticipate any variations of the original physical properties during the tests will be associated with changes in clay content and porosity (e.g., due to clay-washing, stress-induced microcracks, grain removal), and CO 2 -induced clay-swelling.To analyze the effect of these variations on the elastic data collected before and after exposing our samples to CO 2 , we used two approaches: (a) the semi-empirical model proposed by MacBeth (2004) to quantify the weaknesses of sandstones from the pressure sensitivity of the rock frame, and (b) the empirical relationships proposed by Eberhart-Phillips et al. (1989), which enable predicting V P and V S from clay content, porosity and effective stress.By combining these two approaches, we attempt at distinguishing between pure mechanical (e.g., loading/unloading-induced micro-cracking) from chemo-mechanical (i.e., CO 2induced clay-swelling effects) changes in the rock.
In short, MacBeth (2004)'s model assumes the defects of the rock can be described as a set of elements of compliance, which can be evaluated collectively.The pressure (P) dependency of the bulk and shear (K and G) moduli (we use M to refer to any of the moduli, as the same mathematical expression applies to both of them) can be estimated from their background, high-pressure asymptotes M ∞ (K ∞ , G ∞ for their corresponding expressions), as follows: with P M being the characteristic pressure constant (the rock is no longer sensitive to pressure changes beyond this point), and E M a fitting parameter that expresses the sum of cracks and contacts of the unconfined rock.Then, stress sensitivity (S M ) is related to E M , as follows: Journal of Geophysical Research: Solid Earth 10.1029/2024JB029027 FALCON-SUAREZ ET AL.
To complete our assessment, we have also considered the stress sensitivity (S η ) to the P/S velocity ratio coefficient η = (V S /V P ) 2 , expressed as follows: For the modeling, we have assumed an invariable density with increasing pressure, as the voids reduction occurring in consolidated reservoir rocks like the samples used in this study have insignificant impact on the fittings (MacBeth, 2004).

Setting Experimental Conditions
Simulating reservoir stress conditions in the lab is relatively straightforward using triaxial vessels and automatic pressure controllers (e.g., Falcon-Suarez et al., 2020;Hamza et al., 2021;Han et al., 2011;Holt et al., 2000).Temperature control requires more sophisticated setups, which normally include ovens or thermal baths to thermo-control the vessel and fluids (e.g., Munkejord et al., 2020).As most of the sensors are temperaturesensitive, commonly only part of the experimental setup can be subjected to temperature-control, leading to local heat transferences within the setup.These local temperature gradients can play a crucial role when the experiment involves CO 2 , as its properties are highly temperature-dependent in the vicinity of the supercritical point (i.e., scCO 2 is around 31.1°C and 7.38 MPa).Even when we can successfully control the temperature, injecting scCO 2 might damage the equipment.For instance, the high fugacity of the scCO 2 enhances its penetration through rubber materials, like O-rings and sleeves used to isolate the rock from the confining fluid.
The NOC RPL is a 20°C room temperature facility.To find the lab PT conditions that best represent those in the Aurora's reservoir (i.e., P c ∼ 52 MPa, P p ∼ 26 MPa and T ∼ 65°C), we analyzed the PT-gradients for brine and CO 2 bulk modulus (K), density (ρ) and viscosity (η), which are the more relevant properties for elastic waves and electrical data interpretation.Figure 4 shows that K CO2 follows an approximate constant gradient from 26 MPa, 65°C, down to ∼10 MPa, 20°C, with variations of <10% for ρ CO2 and <20% for η CO2 .To assess the deviation of the rock elastic properties at these conditions, we compare the bulk modulus versus saturation curves at the reservoir pore fluid conditions (i.e., P p ∼ 26 MPa and T ∼ 65°C), and at each pore pressure from 5 to 26 MPa (every 1 MPa) for 20°C using Gassmann (GS) and White-Dutta-Odé (WDO) models for the low (in the order of 10s Hz) and ultrasonic (600 kHz) frequency regimes, respectively.For the WDO, we used a CO 2 patchy-sphere of radius 5 × 10 5 m, according to our previous observations (Falcon-Suarez et al., 2016).For the frame and solid particles, we adopted physical and mineral properties of sample ST3 (Table 1).
We found the minimum standard deviation at 10 MPa (SD = 0.04 and 0.05 for GS and WDO curves, respectively).Figure 5 shows how the deviation between both curves lie approximately within our experimental error.Then, assuming an effective stress coefficient n = 1 in the absence of further information, we set the P c = 36 MPa for the CO 2 flow-through tests, to satisfy Aurora site's P eff ∼ 26 MPa.We also assume seawater brine salinity in the reservoir and used a synthetic 3.5% NaCl aqueous solution for the tests.Note that the samples used in this study were extracted from a slightly deeper location than the planned injection point for Aurora.Therefore, we expect no overloading-induced stress microcracks generation on our testing cores when imposing Aurora's stress conditions.

CO 2 Flow-Through (BCFT) Tests
The experimental time of each test varied between 6 and 7 days, resulting in ∼50 pore volume (PV) throughputs for samples ST1 and ST2, and ∼78 PV for ST3.Imbibition stages in ST1 and ST2 lasted shorter due to setup malfunctions, limiting our data interpretation associated with this stage.
The three samples show similar trends for all the geophysical properties (i.e., , and resistivity), but different magnitudes (Figure 6).During the transition from brine to CO 2saturated brine, the three samples evidence small fluctuations in all the properties that lie within the experimental error for each of them.This observation agrees with the bulk resistivity and P-wave velocity results reported by Falcon-Suarez et al. ( 2018) under similar salinity and pressure conditions and slightly higher temperature, which corroborates that our geophysical tools are insensitive to detect dissolved CO 2 fractions in (at least) 3.5% NaCl solutions.
The CO 2 arrival led to a sharp increase in the bulk electrical resistivity above 50% in all the samples, inversely increasing with porosity (i.e., sample ST3 shows the highest variation, followed by ST1).This porosity-resistivity correlation seems to be unrelated to the clay content (i.e., highest for ST1) or mechanical variations, as the very little deformation observed in all the samples (<0.003%) suggest their porosities remained unaltered during the BCFT tests.Similarly, resistivity sharply drops during imbibition, in all cases, fully recovering after 10 PV of brine flow-through based on ST3 data (i.e., the only sample with a complete imbibition record).
V P decreases with CO 2 in the three sandstones, as expected (e.g., Alemu et al., 2013;Kitamura et al., 2014;Lei & Xue, 2009).Samples ST3 and ST1, with lower porosity and higher clay content, show the expected V P trend for sandstones during CO 2 -brine substitution, that is, a sudden V P drop (3%-4%) related to the CO 2 arrival followed by a progressive decrease (down to ∼6%), partially recovered during imbibition (e.g., Alemu et al., 2013;Falcon-Suarez et al., 2018;Kitamura et al., 2014).On the contrary, the velocity of ST2 shows an increase with the CO 2 arrival that is unjustified by the lowering of density alone.This is followed by a progressive decrease as the CO 2 flow-through continues.However, this increase has previously been observed in high porosity sandstones under similar experimental conditions (Falcon-Suarez et al., 2020;Papageorgiou et al., 2018).It has been attributed to the lowering of effective mobility, leading to larger velocity dispersion in the 10%-20% CO 2 saturation range (Falcon-Suarez et al., 2020).V P trends contrast with those observed for V S , which only show the expected increase with CO 2 content for ST3.

Q P
1 and Q s 1 respond equivalently to their respective velocities.The normalized Q P 1 sharply increases with the CO 2 arrival by 65%, 200% and 400% for ST1, ST3 and ST2, respectively.These increments directly correlate with the clay-sized particles fraction of the samples (i.e., ST1 > ST3 > ST2).Note that by normalizing with respect to the Q P 1 of the brine saturated rock (i.e., original condition), we discount the effect of the physical properties of the rock.This observation suggests that clay particles hamper the increase in P-wave attenuation when CO 2 replaces the parental brine in sandstones.Q S 1 decreases with the increasing CO 2 flowing through for ST1 and, particularly, for ST3, while increases for ST2, with absolute values correlating with sample porosity (i.e., ST2 > ST1 > ST3).

Changes in Elastic and Transport Properties After CO 2 Injection
The S-wave particle motion is polarized at 90°to the P-wave particle motion.As a result, V P /V S changes are related to relative strain changes in two orthogonal directions.Therefore, V P /V S ratio is an important parameter that offers information about rock elastic anisotropy, and both the pore fluid composition and pressure.In our samples, V P /V S decreases with increasing P eff in all three samples when saturated with brine (Figure 7), while the opposite occurs when dry (i.e., saturated with compressible air), as previously observed (Wang et al., 2012).
Sample ST-1 shows the most significant V P /V S variation from before to after when dry (Figure 7a), which can be related to microcracks generated during the BCFT tests and oriented perpendicular to the sample axis (i.e., flow direction), as permeability and resistivity trends remain invariable from before to after the BCFT injection test Figure 6.CO 2 flow-through test results in brine-saturated samples ST1, ST2, and ST3.P-and S-wave velocities (V P , V S ), attenuations (Q P 1 , Q S 1 ) and electrical resistivity (ERT) all normalized with respect to their original brine saturated values (subscript 0), together with axial strains.The tests included a long drainage episode (CO 2 replacing original brine in the pores; solid dots), followed by forced imbibition (brine flowed back into the pores partially replacing CO 2 ; empty dots).The ultrasonic properties were measured at a single frequency of 600 kHz (pulse-echo technique), obtained from Fourier analysis of broadband signals.The vertical gridlines indicate the starting point of the CO 2 injection for each sample (ST1, blue; ST2, red; and ST3, black).
(Figure 7b).Sample ST-2 varies both when wet and dry, which can also be related to the development of microcracks.But, in this case, the induced anisotropy seems to be oriented along the sample axis, as permeability (k) increases and the formation factor (F) (e.g., Muñoz-Ibáñez et al., 2019) decreases from before to after CO 2 exposure.In sample ST-3, the trends remain nearly constant after the CO 2 test for all the measured parameters, with a slightly increase in F that indicates a reduction in rock conductivity.Since permeability remains constant, this effect can be related to a potential radial deformation (influencing the radial distribution of electrodes and computation of the electrical resistivity), or a clay-washing effect, as discussed below.

Data Analysis
The imposed stress conditions are below the original stress in the reservoir from where samples ST1, 2 and 3 were extracted.Hence, we should expect minimal loading/unloading stress-induced micro-cracking effects in our samples.Although the degree of damage is mainly conditioned by the original reservoir stress conditions, the majority of the rock defects are developed during coring from wells and linked to (unloading) stress-release effects (Holt et al., 2000).Distinguishing between original and (BCFT) test-induced damages is essential to assess the effect of the CO 2 injection on our samples.
Figure 8 shows that our dry measurements before and after the BCFT tests are well explained by MacBeth ( 2004)'s model.Sample ST1, with the highest clay fraction, exhibits the most significant variation in its (bulk, K, and shear, G) moduli at the target (36 MPa) effective pressure after the CO 2 injection test, followed by ST3.Sample ST3 seems to carry damage developed during the sample extraction/preparation processes, as this is the only sample not reaching asymptotic moduli above 30 MPa pre-CO 2 injection, that is, when most of the cracks and grain contacts forming the so-called soft pore space fraction (MacBeth, 2004) are closed.The sample ST2, with highest porosity and lower clay fraction, accommodates more rapidly the deformation associated with the soft pores, while showing the lowest variation in its frame moduli after the CO 2 injection test.
The stress sensitivity parameters S K , S G , and S η (Table 2) of the three samples agree with the values reported in MacBeth ( 2004) for samples with similar porosities and clay fractions.The sample ST2 shows small variations in all the parameters that suggest little CO 2 -induced damage.For sample ST1, S K and S G decrease post-CO 2 injection, which indicates potential CO 2 -induced microcracking.Also, the drop of S η indicates larger clay fraction post-CO 2 injection (MacBeth, 2004), which agrees with the small decrease observed in porosity (Table 1).In sample ST3, S K and S G decrease, which together with the large increase of P K suggests a preferential closure of the original cracks oriented perpendicularly to wave propagation.Here, again, S η decreases, as well as the porosity of the sample.2.

Table 2
MacBeth ( 2004)'s Model Fitting Parameters and Stress-Sensitivity Factors for the Three Samples (ST1, 2, and 3) The empirical clay-porosity-pressure relationship proposed by Eberhart-Phillips et al. (1989) for saturated sandstones suggests that the sample ST1 post-CO 2 injection would be better explained if considering an increment of the clay-sized fraction (C) (Figure 9), in agreement with our observations above.For sample ST2, the model better fits the observations if a C reduction of 3%-5% pre-CO 2 injection is considered, and could indicate that the original C was partially washed during the brine saturation and initial stages of flowing (i.e., before starting the measurements).This hypothesis agrees with the porosity increase of sample ST2 after the tests.Sample ST3 is poorly explained by this empirical relationship, as C seems to decrease with the increasing P eff .A plausible explanation is that clay-sized particles are predominantly load-bearing in ST3, as opposed to pore-filling only, with clay being more compliant than major grains (e.g., quartz).Hence, the frame porosity would be preferentially affected by the loading, in this case.
The observed changes in the elastic properties of our three samples can be related to CO 2 -induced chemomechanical phenomena.CO 2 can be partially absorbed into clay minerals (Schaef et al., 2014), causing structural changes within the clay layers.Eventually, this reaction would lead to intergranular clay removal and, in turn, to a change (weakening) on the original elastic properties of the rock (Delle Piane & Sarout, 2016).Thus, the apparent increase of C in sample ST1 can be related to swelling-induced CO 2 -absorption, as ST1 contains higher amount of illite and expandable clays (Table 1).In ST3, the dominant clay is kaolinite.Despite being less prone to swelling (e.g., Hamza et al., 2021), kaolinite can absorb CO 2 (Schaef et al., 2014) and lead to volumetric expansion of up to 0.2% at the experimental pore pressure conditions (Zhang et al., 2019).The fact that the porosity only increased in ST2, the sample with lowest clay content, indicates that the CO 2 -altered clay fraction might be washed (even pre-CO 2 injection as previously mentioned), leaving a family of microcracks that potentially closed during the elastic-transport assessment after the CO 2 flow-through experiment.Our results agree with Delle Piane and Sarout (2016)'s observations that the weakening effect is more significant at low P eff , as the clay-related weakening reduces with increasing stress (e.g., Best, 1992).
Resistivity was transformed into the degree of brine saturation (S w = 1 S CO2 ) using the fitting parameters shown in Table 3.The three samples undertook a rapid transition between the fully and partially brine saturation, with very little data from the S w range 1-0.7 during the drainage stage (Figure 10).The velocities during draining are well-explained by the White and Dutta-Odé (WDO) model in samples ST1 and ST3 by assuming patchy saturation bubble sizes of 0.3 and 0.5 mm, respectively; while sample ST2 shows a velocity increase at highintermediate S w during drainage that is poorly explained by the WDO model, for even the extreme bubble size of 3 mm that is close to the experimental wavelength.This local bumping has been previously observed and theoretically predicted (Falcon-Suarez et al., 2020;Papageorgiou et al., 2018).It explains that in high S w partial saturation the effective mobility drops, and velocity dispersion (which increases the velocity) dominates over the fluid effect which lowers the velocity (Falcon-Suarez et al., 2020).In our experiment, this effect is prominent in the highest porosity sample, in agreement with Papageorgiou et al. (2018).
On the contrary, the imbibition stage is adequately explained by the Gassmann-Wood (GW) model in the three cases.During imbibition, brine partially replaces the CO 2 located in the porous medium, leading to a more uniform fluid distribution (e.g., Falcon-Suarez et al., 2020).This fluid homogenization process seems to be less prominent in sample ST3, which might be related to lower pore connectivity in this sample according to its higher formation factor (Figure 7), leading to a slower fluid redistribution.

Modeling Aurora's Site: Time-Lapse Seismic Feasibility Study
Time-lapse (or 4D) seismic feasibility studies aim to look at the sensitivity of the seismic response to changes in subsurface rocks as a result of engineering activities, including fluid extraction and injection activities at targeted formations and over/under-burden formations as a result of geomechanical stresses.The end result is an assessment of whether the subsurface changes can be detected, given the expected changes of the subsurface, the geometry of the experiment, and noise levels.Seismic wave propagation can change in rocks as a result of several factors: changes in pressure, which affects porosity and the rock stiffness, and may even lead to cracking, changes in fluids, and changes in temperature.4D seismic studies primarily involve detecting changes in velocity, or impedance, between baseline (starting conditions) and monitoring (changed conditions) seismic acquisition times, with changes in attenuation not as widely estimated.
CO 2 injection may induce micro-cracking.Seismic wave propagation can be also highly affected, depending on the mechanical and hydraulic properties of the cracks (Song et al., 2020).On the one hand, the mechanical properties affect seismic wave velocity, with the velocity initially increasing with increasing stress due to compaction, up to the yield point when dilation effects start to dominate after which the velocity starts to decrease (Cartwright-Taylor et al., 2022).On the other hand, CO 2 -induced micro-cracking would affect the fluid movement in the porous medium and lead to squirt flow dispersion (Mavko et al., 2009).The dispersion can be evaluated in a very simplistic way using the isotropic squirt model of Mavko and Jizba (1991).This model allows calculating high-frequency saturated V P and V S from dry velocity-pressure data.Mavko and Jizba (1991) express the magnitude of the dispersion at the experimental ultrasonic frequency (where the pore fluid is only partially drained, i.e., stiffer rock frame) as the difference between the experimental velocity and that of the drained rock frame at the low-frequency Gassmann domain.This modeling approach assumes that both Gassmann and unrelaxed (ultrasonic) velocities will be identical at very high P eff , with deviation being linked to the crack-induced energy dispersion.
Plotting together the moduli obtained from Mavko-Jibza and Gassmann models (Figure 11), we observe good agreements between them at high pressure for ST1 and ST2, while the former underestimates the bulk modulus of ST3 by ∼5%.This mismatch can be interpreted as incomplete cracks-closure at the P eff of our experiments for sample ST3.We also observe dispersion at ultrasonic frequency accurately modeled by the WDO patchy saturation model for the three samples, while dispersion at full saturation seems to occur only for sample ST3.
Figure 10 is consistent with squirt flow phenomena in the V P S w relationship.For instance, the increase in velocity in sample ST2 that is characteristic of lowered mobility due to dynamic squirt flow effects more prevalent in intermediate-high saturations.
Keeping this in mind, a good first-order approach to estimate the seismic velocity of the parental reservoir rock formations would be to assume small (5%) and large (20%) squirt flow dispersion between seismic and ultrasonic frequencies and independent of saturation.However, in this case, the tested samples are only representative of the geological depositional environment and depth, so a direct extrapolation of the absolute values could lead to misinterpretation when modeling partial saturations in the reservoir.Instead, to analyze seismic response in Aurora's reservoir during CO 2 injection, and assess the effect of (petro-) physical properties of the rock formation on the results, we have used an alternative approach based on the relative V P difference from highest to lowest brine saturation (i.e., fully saturated, S w = 1, and at the end of drainage stage of the tests, i.e., S CO2,max ).
Our modeling approach combines the available information in Aurora's site (Figure 1) and our experimental data to elaborate two scenarios to assess the effect of rock heterogeneities (clay and porosity) in our interpretation: (a) assuming four different formations within the reservoir (vertically heterogeneous reservoir, Case 1), and (b) considering all the formations as a homogeneous block (vertically homogeneous reservoir, Case 2).For Case 1, we linked the measured properties for ST1, 2 and 3 to the Cook (ST3), Burton (ST1), Johansen (ST2) and Amundsen (ST2) formations at Aurora, based on their average porosity-clay ratio (see yellow frame in Figure 1).Then, we used the experimental percentual change in V P and V S we obtained for ST1, 2, and 3 at the maximum saturation of CO 2 in each case (i.e., S w = 0.62 for ST1, 0.64 for ST2, and 0.56 for ST3) to obtain the sonic (formation average) velocities, while bulk densities at those saturations were obtained from the well information (Figure 1) and fluid (brine and CO 2 ) density at the reservoir conditions (Figure 4).For Case 2, we assumed the four formations considered above behave as ST1, and apply its velocity changes and saturations to all of them.
We conducted a preliminary 4D seismic feasibility study to assess the seismic response to the change in saturation, using the inputs estimated for the heterogeneous and homogeneous cases.We assumed lateral homogeneity and constructed the P-wave reflection coefficient (R) series for zero offset using the following expression: with subscripts i, j denoting interfaces as a function of travel-time for baseline and monitoring conditions.The reflection coefficient series is convolved with a 22 Hz center frequency Ricker wavelet, which is representative of a North Sea seismic experiment (e.g., Mangriotis et al., 2018).
Figure 12 shows the seismic traces at baseline (pre-CO 2 injection) and monitored (post-CO 2 injection) conditions, along with the difference trace (monitored minus baseline).To assess the 4D amplitude response against typical non-repeatability issues, we compute the normalized root-mean-square amplitude (NRMS) defined as follows (Johnston, 2013): We further estimate the theoretical time-shift, calculated as the change in two-way travel time between baseline and monitor conditions, and compare it to time-shift estimated from the cross-correlation method of Rickett et al. (2007) (Figure 12).We perform a further iteration of NRMS calculation after time-shift correction to reflect 4D seismic processing workflows for 4D quantitative interpretation, and a further time-shift estimation after timeshift correction to assess the performance of the time-shift algorithm against a theoretical zero time-shift.
NRMS before and after time-shift application drops from 0.38 to 0.24, for Case 1 (heterogenous), and from 0.74 to 0.42, for Case 2 (homogeneous).Note that typical "good" values of NRMS quoted in the literature range from 0.1 to 0.3 (Johnston, 2013), with further lower NRMS values achievable via seismic permanent reservoir monitoring.We anticipate that further changes due to seismic attenuation will result in a larger 4D difference, suggesting that a 4D signal would be detectable.In terms of estimated time-shifts, 4 ms lies above limits of measured slow-down during gas saturation (MacBeth et al., 2019), suggesting that the time-shift response, similar to the 4D differences, should be detectable above typical noise levels.
The modeling shows a marked difference in the 4D seismic response, both in terms of trace differences and in terms of time-shifts for the two scenarios considered (i.e., Case 1 and Case 2).This observation implies that for quantitative 4D seismic interpretation and, consequently, accurate reservoir dynamic parameter characterization, it is crucial to characterize any reservoir petrophysical heterogeneity where present.

Discussion
Commonly, rock physics models are developed from laboratory observations.Here, we have proceeded in the opposite direction, by using well-accepted rock physics models and theory to set reliable target reservoir PTconditions in our laboratory.Our final goal is to ensure we generate reliable geophysical data during brine-CO 2 fluid substitution tests to properly assess the effect of rock heterogeneities on our geophysical and hydromechanical interpretation of CCS reservoirs-in this particular study, the Aurora's CCS site (a saline aquifer in the northern North Sea).
Uncertainties start when we tried to accommodate the reservoir conditions and needs to our experimental capabilities.Particularly important is the effect of the temperature on the quality of results, as it is a difficult parameter to control in large facilities (such as labs needed for rock physics tests) and might also cause malfunction of some sensors and materials during high-temperature tests.We found that commonly used rock physics models for predicting pore fluid substitution in rocks, such as Gassmann or White-Dutta-Odé models (Mavko et al., 2009), show little difference (i.e., lying within the experimental errors of our geophysical tools; Figure 10) when considering the PT conditions in Aurora's reservoir (P c and P p = 52 and 26 MPa, and 65°C temperature), with respect to those of an equivalent PT point estimated by combining effective stress law and fluid properties to accommodate at our RPL (20°C) temperature.
The main difference between these two (i.e., Aurora and RPL) scenarios is the way the CO 2 is distributed in the pores according to the CO 2 -brine-rock wettability.P p would increase the CO 2 wettability (Al-Khdheeawi et al., 2017;Iglauer, 2017), and we should expect a more homogeneous distribution in the reservoir than in our tests during drainage, but similar to our observations during the imbibition stage.The effect of increasing temperature is less understood, but seems to slightly decrease CO 2 wettability in the presence of quartz and mica (Broseta et al., 2012).Then, in a sandstone aquifer such as Aurora, both the temperature and pressure might counteract one another.In any case, our experimental and modeling approaches covered both drainage and imbibition, providing the high and low boundaries of the elastic properties of the rock for patchy and homogenous CO 2 distribution, which could represent two water-wet extreme scenarios.
The effect of wettability, however, might be negligible with respect to variations in elastic properties due to rock heterogeneities.Our rock samples showed different degrees of rock weakening post-CO 2 testing, which indicates the differences in porosity and clay content highly condition both the hydromechanical response of the rocks and the CO 2 distribution within the pores.According to the high V P /V S ratios shown by the three samples (above 1.8), crack anisotropy seems to be high in all of them (Wang et al., 2012).
We explain the observed changes in the elastic properties of our samples as a combination of different phenomena, including early washing of clays pre-CO 2 injection, chemo-mechanical effects related to CO 2 being partially absorbed onto kaolinite, and an increase of the clay fraction associated with swelling-induced CO 2absorption (Schaef et al., 2014;Zhang et al., 2019).In this regard, we found evidence in our geophysical data of a relationship between the increment in P-wave attenuation when CO 2 arrives in the samples and their clay fractions (Figure 6).In essence, the higher the clay-sized fraction, the lower the Q P 1 variation with CO 2 .This effect can be linked to a reduction of the free CO 2 in the pores due to CO 2 absorption by clay minerals.The quantification of this relationship needs further experiments with samples covering a wider range of clay mineral fractions.Apart from that, the CO 2 -clay induced weakening was geophysically undetectable during the CO 2 flow-through test.
The reasons may include that (a) this effect is masked by the fluid substitution effect; (b) it lies within the measurement error; and (c) it only occurs after washing the sample (i.e., detached clays removed) post CO 2 exposure (when preparing the sample for the elastic-transport assessment after the CO 2 flow-through).
Interpretations derived from the upscaling of ultrasonic waves attributed to the seismic frequencies (used for reservoir monitoring) require careful consideration.On the one hand, the high resolution of the ultrasonic attributes provides insights of the expected changes during fluid substitution associated with textural and mineral heterogeneities.On the other hand, the detection of heterogeneities intrinsically depends on the wave length (λ), which is around five orders of magnitude larger in the field (e.g., Batzle et al., 2006).Basically, variations in rock properties below λ will be overseen and disregarded in the interpretation.However, our results suggest that when CO 2 and brine coexist, even subtle variations in the rock's physical properties can largely affect the distribution of the two pore fluids and the resulting seismic-scale signal.The sensitivity of 4D seismic data to reservoir heterogeneity highlights the importance of laboratory testing of different rock samples representative of the individual formations forming CCS reservoirs.

Conclusions
The experimental and modeling results derived from this study suggest that reservoir rock heterogeneities can largely affect our geophysical interpretation during CO 2 storage operations.
The experimental data collected during CO 2 injection tests in the laboratory show that CO 2 distributes differently in rock samples of the same reservoir formation depending on porosity, clay-sized particles fraction and the clay minerals composition.Changes in P-wave attenuation correlate with clay fraction during the CO 2 injection, together with evidence of CO 2 -induced rock weakening.This provides evidence of fluid-CO 2 -clay coupled chemo-mechanical phenomena during the experiments.
Our 4D proof-of-concept time-lapse seismic modeling shows that disregarding the effect of rock heterogeneities on the partial fluid distribution can lead to significant misinterpretations during CCS operations in terms of both pore fluid distribution and CO 2 quantification.

Figure 1 .
Figure1.Well-log geophysical information from well 31/5-7: P-and S-wave sonic velocities (V P and V S ), bulk density (ρ), porosity (ϕ) and clay-sized particles volume fraction (C) derived from gamma ray (GR) normalized with respect to the considered GR range for each value (i) as C = (GR i GR min )/(GR max GR min ).Solid dots are average values of each parameter within a given geological formation.Yellow frame shows the reservoir section considered for modeling in this study (see main text).

Figure 2 .
Figure 2. Thin sections for the three samples ST1 (a), ST2 (b), and ST3 (c), and results of the three phases segmentation (i.e., grains, pores and clay-sized fractions-estimates are shown next to the segmentation images) on each of them (a.1-c.1) using the machine learning tool Trainable Weka Segmentation, Fiji/ImageJ.

Figure 4 .
Figure 4. Pressure and temperature dependencies for CO 2 modulus (a), density (b) and viscosity (c).Aurora reservoir conditions are marked (square) together with those adopted in the lab (circle).

Figure 5 .
Figure 5.Comparison between the bulk modulus of sample ST3 versus saturation at the PT reservoir (26 MPa, 65°C) and laboratory (10 MPa, 20°C) conditions, using Gassmann (GS) and White-Dutta-Odé (WDO) models.The vertical error bar refers to the estimated experimental error for the ultrasonic measurements and the degree of brine saturation (S w ) derived from electrical resistivity (see main text for further info).

Figure 7 .
Figure 7. Stress dependency of the ultrasonic P-and S-wave velocity ratio, V P /V S , dry and wet (brine saturated), together with absolute permeability (when dry), k, and formation factor (when wet), F, for the three samples (a and b, for ST1; c and d, for ST2; and e and f, for ST3), before (open symbols) and after (solid symbols) being subjected to CO 2 injection.

Figure 8 .
Figure 8. Dry bulk (K) and shear (G) moduli versus effective pressure, data and model fits based on MacBeth (2004)'s model, for the three samples (ST1, 2, and 3) before and after the CO 2 injection tests.See fitting parameters in Table2.

Figure 9 .
Figure 9. Experimental data for wet P-and S-wave velocities before (open symbols, subscript 0) and after (solid symbols, subscript 1) the CO 2 flow-through tests, together with the empirical model fits for saturated sandstones proposed by Eberhart-Phillips et al. (1989).We used the porosity, ϕ, and clay-size volume fraction, C, shown in Table 1, for samples ST1 (a, b), ST2 (c, d), and ST3 (e, f).

Figure 10 .
Figure10.Experimental ultrasonic P-and S-wave velocities (V P and V S , obtained from Fourier analysis of broad band signals at 600 kHz) for drainage (solid dots) and imbibition (empty dots), versus degree of water saturation (S w ), displayed together with the predictions from the White and Dutta-Odé model (WDO; dashed lines) with patchy saturation bubble sizes of 0.3, 3 and 0.5 mm for ST1, 2 and 3, respectively, and Gassmann (GS; solid lines) for 20°C (lab temperature) and 65°C (Aurora's reservoir temperature) scenarios.Error bars displayed every 10 data points.

Figure 11 .
Figure 11.Bulk (K) and shear (G) moduli of samples (a) ST1, (b) ST2 and (c) ST3 estimated from ultrasonic data at fully saturated conditions in brine, together with the corresponding Gassmann and Mavko-Jizba (squirt) models.

Figure 12 .
Figure 12. 4D seismic modeling.(a) Seismic amplitude baseline and monitored traces (B and M, referring to pre-and post-CO 2 injection conditions; see main text), (b) M-B difference and (c) the time shift; (d, e, f) same plots after time-shift correction.Solid and dashed lines show results for the Case 1 and Case 2 (i.e., vertically heterogeneous reservoir and vertically homogeneous reservoir), respectively (see main text).

Table 1
Physical Properties and Mineralogy of the Core Samples Used in the Tests

Table 3
(Juhász, 1981) Parameters Used for the Rock Physics Modeling, and Fitting Parameters for the Transformation of Resistivity Into Degree of Saturation Using the Modified Archie's Relationship to Account for the Clay Effect(Juhász, 1981) FALCON-SUAREZ ET AL.