Drill core processing

HQ core (63.5 mm diameter) samples of the caprock interval, recovered from well CO2W55 near the town of Green River, Utah, were slabbed and then sliced parallel to bedding at ∼3–5 mm resolution, using a steel blade rock saw. Approximately 2–5 g of air-dried sample was powdered to ≤100 μm size in a ball mill using an agate container and balls, and subsamples were taken for individual analyses. Samples for XRD measurements were crushed and milled separately. A comparative sample profile from the base of the Carmel Formation, in a region where CO 2 is absent, was collected from sections of core BH2 recovered during drilling of the Bighole fault44, located ∼35 km north east (NE) of Green River. This comparison profile was prepared and analysed for quantitative mineralogy and carbonate O- and C-isotopes using the same methods as for the CO2W55 core interval (Supplementary Fig. 3).

XRD measurements

XRD measurements to determine mineralogy (Supplementary Table 1) have been performed at RWTH Aachen University as stated in ref. 45, except the counting time is 20 s for each step of 0.02° 2θ recorded from 2° to 92° 2θ. Rock samples are crushed manually in a mortar with care taken to avoid strain damage and crushed material together with an internal standard (Corundum, 20 wt%) is milled in ethanol with a McCrone Micronising mill (15 min). Quantitative phase analysis is performed by Rietveld refinement using BGMN software, with customized clay mineral structure models46. The precision of these measurements, from repetitions, is better than 0.1 wt% for phases of which the content is above 2%. The accuracy cannot be determined because of the lack of pure clay mineral standards, but is estimated to be better than 10% (relative). Mineral compositions relate to the crystalline content of the analysed samples.

Stable isotopes

Carbonate mineral δ13C and δ18O were determined at the University of Cambridge on duplicate subsamples of 300–500 μg using a Thermo Gas Bench attached to a Thermo MAT 253 mass spectrometer in continuous flow mode, with an analytical precision of ±0.12‰ and ±0.20‰, respectively, (2σ). δ13C and δ18O results are reported, relative to the VPDP and VSMOW standards, as averages of these duplicate analyses, with error bars that are the 2σ s.e. of duplicate averages (Supplementary Table 2 and Supplementary Fig. 1).

Sr-isotopes

87Sr/86Sr of rock leachates and residues were determined at the University of Cambridge following ref. 47 (Supplementary Table 2 and Supplementary Fig. 1). The internal standard NBS 987 gave 0.710263±0.000009 (1σ) on 158 separate measurements made during the course of these analyses. Strontium blanks were always <250 pg and negligible for the Sr concentration of these samples.

EMPA and scanning electron microscopy analyses

EMPA of mineral compositions were determined using a Cameca SX-100 electron microprobe, using energy-dispersive spectrometry at the University of Cambridge (15 kV, 10 nA; beam diameter 5 μm) with fayalite, rutile, corundum, periclase and pure Co, Ni, Mn, Cr, Zn and Cu standards (Supplementary Table 4). Spectra were collected with a PGT prism 2,000 ED detector and the data reduced with the PGT excalibur software. Back-scatter electron and secondary electron images were obtained using the JEOL 820 scanning electron microscope in the Department of Earth Sciences, University of Cambridge, with an accelerating voltage of 20 kV and 1 nA beam current.

Porosity measurements by N 2 physisorption

The surface area and pore-size distribution of sample NPS-069 was investigated using nitrogen gas adsorption techniques at RWTH Aachen University as stated in ref. 35 (Supplementary Table 3 and Supplementary Fig. 6). Specific surface area and total pore volume were determined from nitrogen gas adsorption at 77.3 K, by means of the static-volumetric method, using a Micromeritics Gemini VII 2390t. Samples were prepared by gentle manual crushing of the supplied core sample, with minimal use of energy. The crushed material was manually sieved, the measurements were performed on the 63–400 μm size fraction.

About 1 g of crushed sample material was outgassed in a Micromeritics VacPrep 061 for 12 h at room temperature and then heated to 130 °C, under vacuum, for at least 12 h. Adsorption was measured at 42 relative pressure steps between 0.01 and 0.995, and desorption at 29 relative pressure points between 0.995 and 0.1. The nitrogen saturation pressure (P 0 ) was determined separately for each relative pressure point. Sorbate–sorbent equilibrium was assumed when the pressure change over a 10-s interval was <0.01% of the average pressure during the latter interval. The multipoint Brunauer–Emmett–Teller theory (BET) method was applied to quantify the surface area of the analysed sample. A cross-sectional area of the nitrogen molecule of 0.162 nm2 was used (ISO 9277:2010). Total pore volume was determined by means of the Gurvich rule at an interpolated relative pressure of 0.995, assuming the density of adsorbed N 2 equals that of liquid N 2 at the boiling point (28.8 mol l−1). Pore volume and area distribution were calculated by means of the Barrett–Joyner–Halenda method in accordance with the DIN 66,134 norm. The Harkins–Jura equations with Faas correction were used for calculation of the statistical thickness curves. The reported Barrett–Joyner–Halenda data are calculated from the adsorption branch of the isotherms, assuming half of the pores are open at both ends. Although desorption isotherms are generally advised for pore size distribution (PSD) calculations, these do not give useful data for shales and other geological materials because of the very strong overprint by the tensile strength effect of N 2 at ∼40 Å. Repeated measurements on standard materials demonstrated that the accuracy and precision of the methods described above is better than 2% (2σ).

Small angle neutron scattering

Sample porosity and specific surface area was analysed using SANS and VSANS techniques (Fig. 3, Supplementary Figs 4 and 6, and Supplementary Table 3). Experiments were carried out using the instrument KWS-1 (SANS) and KWS-3 (VSANS) operated by the Jülich Center for Neutron Science at Heinz-Meier-Leibnitz Zentrum in Garching, Germany. Carmel caprock samples were cut parallel to bedding, fixed on quartz glass and polished to a thickness of 200 μm for measurements. Samples were dried at room temperature and measurements performed under ambient pressure and temperature conditions. The target area on the sample was defined by a cadmium mask with an 8 mm diameter window.

For (V)SANS measurements, a collimated neutron beam is elastically scattered by the sample48,49. Position-sensitive detectors measure the scattering intensity I(Q) as a function of the scattering angle, which is defined as the angular deviation from the incident beam. The momentum transfer Q (nm−1) is related to the scattering angle θ by Q=(4π/λ)sin(θ/2), where λ is the wavelength of the neutrons. Thus, the size range of features accessible with neutron scattering depends on the neutron wavelength λ and the range in the scattering angle θ. SANS data at KWS-1 were collected at wavelengths of λ=0.69 nm with a wavelength distribution of the velocity selector Δλ/λ=0.10 (full width at half-maximum). Measurements were performed at sample-to-detector distances of 19.7, 7.7 and 1.7 m, covering a wide Q-range of 0.02–2.6 nm−1. The detector was a 6Li glass scintillation detector with an active area of 60 × 60 cm2. VSANS data at KWS-3 were collected at λ=1.28 nm, Δλ/λ=0.2 and a sample-to-detector distance of 9.5 m, covering a Q-range from 0.024 to 0.0016, nm−1. As for KWS-1, a 6Li scintillation detector was used but with a detector diameter of 9 cm. Hence, pore radii for the combined SANS and VSANS measurements range between 1 nm and 1.5 μm (r≈2.5/Q). Instrument data analysis and background subtraction was carried out using the QtiKWS software provided by the Jülich Center for Neutron Science http://iffwww.iff.kfa-juelich.de/pipich/dokuwiki/doku.php/qtikws. During background subtraction, the lower pore sizes were cut off at 12 Å, to remove artefacts arising from Bragg scattering from ordered stacking of clay minerals. For the combined SANS and VSANS measurements, this results in scattering intensity I(Q) (nm−1) versus the momentum transfer, Q (nm−1) relationships on nine samples (Supplementary Fig. 4).

The integral of the scattered intensity in reciprocal space in a two-phase system is the Porod invariant Q inv from which the porosity ϕ is obtained directly48,50:

We assumed that the scattering intensity of pore features is directly proportional to the scattering contrast between matrix and pore48:

Here, and are the coherent scattering length densities (SLDs) for neutrons for the two phases, shale matrix and air (pore), respectively. The terms P(Q) and S(Q) denote the form and structure factors, for which analytical expressions exist for different geometries of scatterers, including mass and surface fractals48. The SLD Δp* of each mineral was calculated as:

where b i and M i are scattering length and atomic mass of the ith element in the mineral, d i is the mass density of the ith mineral and N A is the Avogadro number. The terms ϕ and V p are the volume fraction of the dispersed phase and the volume per scatterer, respectively. Scattering length densities were calculated using the U.S. National Institute of Science and technology (NSIT) SLD calculator (http://www.ncnr.nist.gov/resources/activation/). As the scattering contrast between the shale matrix and the pores is large, all scattering is attributed to nano-pore features. SLD values for the shale matrix studied range between 3.2 and 3.4 × 10−14 m−2; for air, it can be considered to be zero.

The specific surface area of surface fractals scales with the length scale r as48,51:

where ρ is mass density, D s is the surface fractal dimension and

Experimentally determined I(Q) curves were modelled, after background subtraction, using the polydisperse hard sphere model implemented in the software code PRINSAS51. Porosity, pore volume and specific surface area values obtained from this analysis are summarized in Supplementary Table 2; intensity, I(Q) and pore frequency distribution f(r) plots for all nine samples analysed are shown in Supplementary Fig. 4. Effective diffusion values D e (m2 s−1) were calculated from the fractal dimensions derived from the (V)SANS data (Fig. 3c). Typically, D e is related to the diffusion coefficient in water, D w and the diffusion accessible porosity, ϕ, by

where τ is the tortuosity. Using fractal dimensions derived from (V)SANS measurements; however, the effective diffusion coefficient is related to the diffusion accessible porosity by an analogue empirical power law formulation of the form:

where m is an empirical exponent. In natural porous media the power law form (equation 13) is indicative of the fractal geometry of the pore–solid interface. These fractal dimensions of the pore surface and pore volume can be quantified from (V)SANS analysis. In such a sinuous capillary bundle model, the sinuous nature of the bundle reflects the roughness of the pore–solid interface, with fractal properties within some scale limits. The power law form (equation (13)) can thus be related to the porosity of a volume element of size l 2 , between the lower and upper limits of the fractal region, l 1 and l 2 , and to the pore volume fractal D v , as

where, for embedded dimensions of three52, 2≤D v ≤3. The tortuosity(Fig. 3c) of the fractal bundle in a volume element is related to the fractal dimensions of the pore–solid surface, D s , and pore volume, D v , as

Substituting equation (15) into equation (12) gives

where the fractal dimensions of the pore surface and pore volume can be quantified from (V)SANS analysis as the slope of Q versus I(Q) and r versus f(r), respectively, in the Q-range 10−4 Å−1≤Q≤10−2 Å−1.

Hydraulic conductivity measurements

A section of core from a red siltstone unit of the lower Carmel Formation (depth: 180.4–180.7 m) was selected during drilling and preserved under mechanical confinement on-site for later analysis in the laboratory. Intrinsic permeability of subsections of the core interval were measured in a bespoke permeameter designed and built by the British Geological Survey. Each core plug sample was carefully manufactured on a machine lathe to minimize damage and desaturation, and to produce samples of known dimension. Samples were then placed inside a single closure pressure vessel and subject to in situ stress and porewater pressure conditions, to hydrate the sample before testing. Permeability was measured by the application of a fixed pressure gradient across each core, while simultaneously measuring flux in and out of the sample. Volumetric flow rates into and out of the core were controlled or monitored using a pair of ISCO-260, Series D, syringe pumps operated from a single digital control unit. The position of each pump piston is determined by an optically encoded disc graduated in segments equivalent to a change in volume of 16.6 nl. Movement of the pump piston is controlled by a micro-processor that continuously monitors and adjusts the rate of rotation of the encoded disc using a DC motor connected to the piston assembly via a geared worm drive. This allows each pump to operate in either constant pressure or constant flow modes. A programme written in LabVIEW elicits data from the pump at pre-set time intervals. Testing was performed in an air-conditioned laboratory at a nominal temperature of 20 °C.

Darcy’s law gives the following relationship for the conductivity, K,

where Q is the steady-state flow (m3 s−1), p w is the density of water (kg m−3), g is the acceleration due to gravity (m s−2), A s is the cross-sectional area of the test sample normal to flow (m2), ΔL s is the sample length (m) and ΔP is the pressure drop along the sample (Pa). The equivalent permeability term (k) is given by

where μ w is the viscosity of water (1.002 × 10−3 Pa s). Measurements on the four selected core intervals gave permeabilities between 5 × 10−22 and 1.2 × 10−21 m2 with an average value of 9.1 × 10−22 m2 (1σ=2.3 × 10−22 m2).

Reaction rates for dolomite and haematite

Dolomite dissolution rates are given by Pokrovsky et al.53 who modelled experimental determinations of dolomite dissolution rates by various pH-dependent surface processes. At the intermediate pH conditions and CO 2 partial pressures in the cap rock, the predicted reaction rates are between 5 × 10−7 and 10−5 mol m−2 s−1. dos Santos Afonso and Stumm41 measured haematite dissolution rates as a function of pH and H 2 S concentrations. For the pH and estimated partial pressures of H 2 S in the caprock, the predicted dissolution rates lie between 10−10 and 10−8 mol m−2 s−1.

Reactive transport modelling with analytical solution

Advective flow rates, (ω 0 ϕ) where ω 0 is fluid velocity and ϕ porosity, through the caprock are estimated to be <6 × 10−12 m s−1 given permeabilities of ∼10−21 m2 and a piezometric pressure gradient of ∼104 to 4 × 106 Pa m−1 taking the hydraulic head in the Navajo54 either across the whole Carmel Formation or the 16 cm basal seal. With effective diffusivities (D e ) of components in the fluid phase of ∼10−11 to 10−12 m2, the Péclet number (ω 0 ϕh/D e ) for advective–diffusive transport of a compatible component would be in the range 10−7 to 10−2, values that imply advective transport is negligible consistent with the analysis of the caprock transport. We therefore model reactive transport in the caprock driven only by diffusion.

The equation describing one-dimensional diffusive transport with mineral reaction of solute concentration C (mol m−3) with respect to distance x (m), time t (s), effective diffusivity D (m2 s−1), porosity ϕ, tortuosity, τ2, mineral reaction rate k f (m s−1) and surface area α (m2 m−3) can be written as39

where (C eq −C) describes a linear dependence of the reaction rate to the difference between the solute component concentration, C (mol m−3) and its concentration at equilibrium, C eq , with the reacting mineral. The mineral molar concentration in the rock, ϕ s (mol m−3), then varies according to

Making the transformation to dimensionless variable by

where l (m) is a transport distance taken as the displacement of the reaction front, D e is the effective diffusion coefficient ϕD/τ2 and C 0 is solute concentration in the input fluid. Substituting equation (21) into equations (19) and (20) gives

dependent on the dimensionless constant, the Damköhler number N D , defined by

For reaction rates of dolomite (k R ) between 10−6 and 10−8 mol m−2 s−1 (note k f =k R V S where V S is the molar volume of dolomite) and diffusion coefficients in the range 1 × 10−11 to 5 × 10−12 m2 s−1, Damköhler numbers are in the range 5 × 107 to 9 × 109 and fluids will be effectively in local equilibrium with minerals. For haematite with reaction rates in the range 10−8 to 10−12 mol m−2 s−1, Damköhler numbers are in the range 2 × 103 to 3 × 107 and haematite will be close to equilibrium wherever it is in contact with CO 2 transported by diffusion.

On initiation of diffusion, dissolution of the reactive mineral takes place at a rate, which decreases away from the inlet such that the mineral is exhausted after time τ 0 (s) given by

or in dimensionless time,

At times greater than τ 0 , the reactive mineral is exhausted over distances x≤l. Where concentrations in the fluid are small compared with those in the solid and porosities are small, the approximate quasi-stationary state solutions to equations (19) and (20) are valid39. Solutions for fluid compositions and solid mineral modes, assuming mineral surface area remains constant, following ref. 39, in dimensionless coordinates for τ 0 >0, are (note l′=1)

where ΔC 0 =C eq −C 0 . ΔC 0 for Fe 2 O 3 in reaction (1), between the inlet solution and a solution at equilibrium with haematite, has been calculated using PHREEQC42 as 2.95 mol m−3 for an inlet solution with the average composition of the fluids sampled downhole in the Navajo Sandstone, with a range from 1.9 to 3.95 mol m−3 reflecting the range in pH of 5.3–5.1 of the downhole fluids23. It is noteworthy that the Damköhler number, N D =(ql)2, with q (m−1) defined by Lichtner39 as the exponential constant giving the length scale over which fluid returns to equilibrium with the initial mineralogy. The time for the reaction front to migrate to distance, l, (l′=1) is given by solution of

or in terms of dimensional constants

Multicomponent reactive transport modelling

PHREEQC modelling was conducted using the llnl.dat thermodynamic data base thermo.com.V8.R6.230 (ref. 42). A 15-cm-long one-dimensional reactive–diffusive model comprising 30 cells of 5 mm length was constructed. The initial mineralogy (mol l−1) was calculated from quantitative mineralogy determined by XRD, using pore volumes measured from SANS of the unaltered portion of the caprock. The invading pore fluid chemistry was based on analyses of the reservoir fluids (Supplementary Table 4) and the initial caprock pore fluid chemistry was taken to be a fluid in equilibrium with the caprock mineralogy, with pCO 2 , pO 2 and salinity estimates for typical Jurassic marine shales. The initial redox state of the invading fluid was defined using the SO 4 2−/H 2 S redox couple. Models were run assuming local fluid–mineral equilibrium. A constant D e value of 5 × 10−12 m2 s−1 was used for all aqueous species. The model timescale was 125,000 years with a time step of 7 days. The thermodynamic database was supplemented with data for Fe-dolomite (CaMg 0.9 Fe 0.1 (CO 3 ) 2 ) calculated using ideal mixing laws and the thermodynamic data for ankerite taken from ToughReact TherAkin8.dat.

A second model was constructed comprising two 25,000 year pulses of CO 2 -charged brine injection separated by a 75,000 year period where a hypothetical CO 2 -poor brine was diffused into the model cells (Fig. 4b). The durations of each pulse were chosen to mimic the timescale of periodic degassing of CO 2 previously document for the local fault systems26.

Data availability

The authors declare that the data supporting the findings of this study are available within the article and its Supplementary Information files.