Observational evidence confirms modelling of the long-term integrity of CO2-reservoir caprocks

Storage of anthropogenic CO2 in geological formations relies on a caprock as the primary seal preventing buoyant super-critical CO2 escaping. Although natural CO2 reservoirs demonstrate that CO2 may be stored safely for millions of years, uncertainty remains in predicting how caprocks will react with CO2-bearing brines. This uncertainty poses a significant challenge to the risk assessment of geological carbon storage. Here we describe mineral reaction fronts in a CO2 reservoir-caprock system exposed to CO2 over a timescale comparable with that needed for geological carbon storage. The propagation of the reaction front is retarded by redox-sensitive mineral dissolution reactions and carbonate precipitation, which reduces its penetration into the caprock to ∼7 cm in ∼105 years. This distance is an order-of-magnitude smaller than previous predictions. The results attest to the significance of transport-limited reactions to the long-term integrity of sealing behaviour in caprocks exposed to CO2.

C arbon capture and storage will form an essential part of the technologies needed to reduce anthropogenic CO 2 emissions if costly climate change is to be avoided 1 . CO 2 , separated at power stations and industrial plants, is compressed and injected into saline geological reservoirs as a supercritical fluid 2 . The CO 2 is less dense than the saline brines at typical formation temperatures, driving buoyant migration, and will need to be retained by impermeable caprocks as the primary mechanism for ensuring effective storage. Subsidiary processes including dissolution of CO 2 in formation brines, trapping of CO 2 by capillary forces and precipitation of carbonate minerals will aid storage security over the 10 4 year timescales needed to avoid climate impacts 3 , but operation of these processes is dependent on the initial retention under the caprock. The caprocks in sedimentary successions suitable for CO 2 storage will most commonly comprise clay-rich mudstones or shales that act as barriers to CO 2 migration by two mechanisms 4 . Their low permeabilities (o10 À 19 m 2 ) restrict fluid fluxes to very low rates. Further, capillary entry pressures for CO 2 of 0.5-5 MPa 4 will restrict penetration of CO 2 , allowing transport of CO 2 only by sluggish diffusion in the brine phase, which will be further retarded by the complex geometry of the pore networks. The addition of CO 2 forms acid brines that react with clay-rich caprocks to drive dissolution and precipitation of minerals [5][6][7] , with consequences for the evolution of porosity and the geomechanical integrity of the seals 8 , within regions penetrated by the diffusing CO 2 . The initial dissolution of rapidly reacting carbonate and oxy-hydroxide phases will buffer the pH of CO 2 -saturated brines to B5 followed by the more sluggish dissolution of silicate and phyllosilicate minerals, which further increase pH and cause re-precipitation of carbonate phases 5,9,10 . However, it is uncertain as to whether these fluid-rock reactions may lead to mineral precipitation-induced self-sealing phenomena that limit the pervasion of the CO 2 , as predicted by some numerical simulations 5,[9][10][11][12][13] , or self-enhancing mineral dissolution and porosity generation, which generate a continuous increase in caprock transport properties, as observed in some laboratory experiments [14][15][16][17] . Direct observations of caprock materials exposed to CO 2 over the timescales requisite for storage are critical to address this uncertainty, where the coupled reactive-transport phenomena are integrated over the appropriate temporal and spatial scales.
Fundamental uncertainties in the prediction of coupled reactive transport in low permeability clay-rich caprocks include a poor understanding of the relationship between reaction-induced changes in porosity and the pore-network structure, and the consequences for the rates of solute transport 18 , the uncertain mineral surface areas and kinetics of the mineral reactions in natural settings 19 , the rates of which may be limited by the intrinsic surface reaction rates of the minerals or by rates of solute transport 20 , and the reaction pathways in natural settings, including the role of acid redox reactions as a CO 2 sink 21 , which will retard the rates of diffusive CO 2 transport.
We address this fundamental gap in our knowledge of the long-term impacts of CO 2 -charged fluids on caprock integrity by examining caprock recovered from a natural CO 2 reservoir. Although this reservoir is at much shallower depths and lower pressures than needed for geological carbon storage, the reactions between minerals and CO 2 -saturated brines are little affected by pressure. Mineralogical and petrophysical measurements are used to assess the impacts of CO 2 -promoted reactions on the caprock pore-network transport properties. The measurements reveal that the CO 2 -charged mildly reducing brines only removed haematite from the basal 7 cm of the haematite-bearing claystone cap rock. This data is used to constrain analytical and numerical reaction-diffusion models, which confirm that the fluid-rock reactions retard CO 2 transport by an order of magnitude, and that the timescale of alteration is B10 5 years, comparable with geochronological constraints on the duration of CO 2 storage in the reservoir. The small but transient increases in porosity generated by the mineral reactions do not significantly enhance reaction front velocities. The results attest to the significance of transport-limited reactions to the long-term integrity of sealing behaviour in caprocks exposed to CO 2 .

Results
Mineral reactions and petrophysical changes in the caprock. Core and fluid samples were collected from a sequence of CO 2 -charged reservoirs and intervening caprocks during scientific drilling within the footwall of the Little Grand Wash Fault, Green River, Utah 22 (Fig. 1). The 325-m deep drillhole transected reservoirs of CO 2 -charged brine in the Middle Jurassic Entrada and Lower Jurassic Navajo Sandstone, and the intervening Middle Jurassic Carmel Formation caprock. The geochemistry of fluid samples, collected at formation pressures through the Navajo Sandstone reservoir, tracks contemporary filling via artesian flow of CO 2 and CO 2 -saturated brines through the faults 23 . The reservoir fluids are mixtures of Na-Cl-SO 4 brine and meteoric groundwater with low pH (5.1-5.3), reducing (0 to À 50 mV), SO 4 -rich and contain trace quantities of H 2 S and CH 4 . The migrating CO 2 and CO 2 -saturated brine are thought to originate from an accumulation of supercritical CO 2 within Carboniferous strata at depth 24 . U-Th geochronology of carbonate veining in associated faults attests to CO 2 out-gassing regionally for ca. 400 ka and locally for at least 114 ka, providing an independent constraint of the duration over which the Carmel Formation caprock has been exposed to CO 2 -charged brines 25 . Cycles in the chemistry of the fault-hosted carbonate veins have been attributed to cyclic charging of the shallow reservoirs with multiple pulses of CO 2 -rich fluids, which coincide with periods of crustal unloading during regional deglaciations 26 .
The Carmel Formation constitutes a regional seal between the Navajo and Entrada aquifers. It comprises a 50-m-thick complex package consisting of three major lithofacies; interbedded, unfossiliferous red and grey shale and bedded gypsum, red and grey claystone/siltstone and fine-grained sandstone 23 . These are interpreted as marine sediments deposited in quiet, subtidal conditions under the influence of periodic hypersaline water 27 .
In the Green River drill hole, a basal 16-cm-thick claystone of the Carmel Formation acts as a caprock to the CO 2 -rich brines in the underlying Navajo Sandstone. This claystone was subsampled at B4 mm vertical resolution and analysed for quantitative mineralogy, carbonate d 18 O and d 13 C isotopic compositions, 87 Sr/ 86 Sr ratios, bulk mineral surface area and porosity (Fig. 2a, Methods, Supplementary Fig. 1 and Supplementary Tables 1-3).
The primary mineralogy of the unaltered portion of the red claystone comprises illite, quartz, dolomite, K-feldspar and haematite ( Supplementary Fig. 2), with 450% clay-sized particles. This basal claystone is typical of oxidized redbed shallow marine and continental mudstones, forming a representative analogue to other important reservoir seals such as the North Sea Mercia Mudstone 28 and the Formations in the Triassic Keuper 29 .
The basal 7 cm of the originally red claystone, at the caprockreservoir interface, is distinctively bleached white-yellow, reflecting quantitative removal of haematite (Fig. 2a), with loss of the primary dolomite and precipitation of ankerite-dolomite, sulphide minerals and gypsum ( Fig. 2a and Supplementary Fig. 1A). Carbonate, sulphate, sulphide and clay mineral compositions have been analysed using Electron microprobe analyses (EMPA; Methods and Supplementary Table 4). Large-scale bleaching of exhumed redbed sandstones and siltstones in the Green River location has previously been related to flow of reducing CO 2 -bearing fluids [30][31][32] . A comparable basal unit from the Carmel Formation from a drill hole free from CO 2 33 km north west (NW) of Green River shows no discernible change in mineralogy or O-and C-isotope ratios as a result of reactions at the base of the caprock ( Supplementary Fig. 3).
The observed petrological changes are consistent with a series of interdependent reactions involving CO 2 -and H 2 S-promoting dissolution of haematite following the stoichiometry Dolomite dissolution was coupled with the growth of ankeritedolomite, with Fe 2 þ derived from the acid-reductive dissolution of haematite ( Fig. 2a and Supplementary Fig. 2A Pore network structures at 4 mm intervals along the reaction profile were characterized using N 2 -BET and small and very small angle neutron scattering techniques ((V)SANS; see Methods). The porosity profile shows a remarkable increase in porosity directly upstream of the haematite dissolution front, but further upstream precipitation of ankerite-dolomite, pyrite, gypsum and illite cause a progressive decrease in porosity to values below that of unreacted caprock ( Fig. 2b and Supplementary Table 3).
The tortuosity (t 2 ) of the pore network at each sample point in the profile (Fig. 2c) was calculated from the surface and volume fractal dimensions of the pore structure, measured by the (V)SANS data (Supplementary Fig. 5 and Methods).
This allows the calculation of the effective diffusivity of a species, D e (m 2 s À 1 ) in the tortuous pore network, which is related to the diffusion coefficient in water, D w , and the diffusion accessible porosity, f, by The molecular diffusivity of CO 2 and H 2 S in water at 20°C is B2 Â 10 À 9 m 2 s À 1 (ref. 33). The calculated porosities, surface fractal dimensions, D s , and pore volume fractal values, D v , (Supplementary Table 3) imply a t 2 of B37 and effective diffusivities for CO 2 and H 2 S of 5 Â 10 À 12 m 2 s À 1 in the unaltered portions of the caprock (Fig. 2c). The calculated porosities agree well with the porosity distribution measured by N 2 -BET over the same size range for sample NPS-069 ( Supplementary Fig. 6). Clarkson et al. 34 and Bertier et al. 35 discuss comparison of porosities calculated from (V)SANS data with N 2 physiosorption data for comparable shales and emphasize the limitations of the physiosorption data. The calculated diffusivities are at the lower limits of experimentally determined measurements for tortuosity and CO 2 effective diffusivities in shales 18,[36][37][38] , and are consistent with the relative high clay content of these samples. It should be noted that our calculations assume that the porosity is interconnected for pores over the length scales sampled (1.2-1,500 nm).
In the altered portion of the caprock the estimated tortuosity decreases to 7 and effective diffusivities increase to 3 Â 10 À 11 m 2 s À 1 , because pore connectivity increases and the roughness of the pore-solid interface decreases due to mineral dissolution (Fig. 2c). The diffusivities decrease to 1.9 Â 10 À 11 m 2 s À 1 further upstream following mineral reprecipitation, but the tortuosity of the pore network does not significantly increase, suggesting that there is a non-recoverable impact of mineral dissolution on rock transport properties.
Modelling reactive transport in the caprock. The mineralogical profiles record the propagation of mineral-fluid reaction fronts within the caprock lithology, driven by the upward diffusion of CO 2 and H 2 S. Modelling of the reactive transport enables examination of the key processes governing the reactive diffusion of CO 2 including the timescale and the ratio of transport rates to Mancos shales mineral reaction rates. Upward advective transport by fluids must be negligible given the low permeability of the unreacted caprock of B10 À 22 m 2 (Methods), which equates to a hydraulic conductivity of B10 À 15 m s À 1 , and thus a Péclet number of B10 À 3 m, over the length scale of the alteration. An analytical solution to one-component diffusive transport with mineral dissolution rate described by linear kinetics 39 captures the important physics and informs the essential approximations necessary for reactive transport modelling with multiple components. On initiation of diffusion, the reactant mineral volume fraction decreases until the phase is exhausted at the base of the caprock. From this time, t 0 , a reaction front migrates downstream away from the base, where t 0 is given by, mineral surface area, V s is the mineral molar volume (m 3 mol À 1 ) and f 1 s is the volume fraction of the reactant mineral phase initially present. DC 0 ¼ C eq À C 0 where C eq is the solute concentration (mol m À 3 ) in equilibrium with the reactant mineral phase and C 0 is the solute concentration in the infiltrating fluid (for full derivation of equations after ref. 39, see Methods), which can be modelled in terms of consumption of CO 2 , H 2 S or Fe 2 O 3 , for which the solutions are identical and fixed by the stoichiometry of reaction (1). The position of the reaction front, l (m) at time t (s) is given by solution of where q (m À 1 ) is given by where D e is the effective diffusion coefficient defined in equation (2). At t4t 0 , the variation of reactant mineral volume with time and distance (x, m) downstream of the reaction front is given by In the absence of advection, solutions for the position and velocity of the reaction front (equation 4) are divided into two regimes and the broadening of the front depends dominantly on the kinetics of the reaction 40 . If ql41 (that is, the Damköhler number (ql) 2 is large), reaction rates are fast compared with diffusion rates and the propagation rate depends on the diffusivity of the species driving the reaction and the reaction stoichiometry. For qlo1, the position of the reaction front is a linear function of time and the velocity of the front depends on the reaction rate constant and mineral surface area. Figure 3 illustrates such solutions, calculated for transport of Fe in the fluid phase, for time, t, as a function of haematite dissolution rate, k R (mol m À 2 s À 1 ) given the reaction stoichiometry of equation (1) and contoured for the effective diffusivity of the caprock, D e , estimated from pore volume and fractal pore network modelling of SANS data.
The haematite mode profile (Fig. 2a) primarily reflects dissolution described by equation (1), but potentially with additional precipitation of haematite downstream of the reaction front as a result of multicomponent transport ( Fig. 2a and c.f. ref. 30). The geometry of the mineral dissolution front, characterized by the spatial variation in haematite modes, is described by equation (6). The curvature of this profile, for an assumed constant effective diffusivity at and upstream of the dissolution front, depends only on the kinetics of the reaction. A least-squares best fit gives a value for the exponential constant q in equation (4) of 228±31 m À 1 (1s). If the profile is augmented by haematite precipitation, q might be as small as B37 m À 1 . Given that l ¼ 0.07 m (Fig. 2a), ql is 42 and reaction rates are fast compared with diffusive transport. For D e ¼ 5 Â 10 À 12 m 2 s À 1 (Fig. 2c), the reaction stoichiometry of equation (1), a haematitespecific surface area calculated from a volume-weighted total SANS surface area and inlet fluid compositions from downhole sampling 23 (Methods), these values of q imply a timescale between 3,000 and 10,000 years, and haematite dissolution rates in the range 1 Â 10 À 10 to 5 Â 10 À 9 m 2 s À 1 (Fig. 3). The latter compare well with laboratory estimates for haematite dissolution rates under acidic conditions in the presence of low concentrations of dissolved H 2 S 41 (Methods).
The modelling above only considers transport of Fe 2 O 3 and dissolution of haematite as the major component buffering fluid pH and oxidation state. The mineralogical profiles ( Fig. 2a and Supplementary Fig. 1A) show that dissolution/precipitation of dolomite, ankerite-dolomite, sulphate, sulphide, K-feldspar and illite also occur. The reactive diffusion was modelled with multiple components and phases by numerical simulation using PHREEQC 42 (Fig. 4, Methods). A 15-cm-long one-dimensional reactive-diffusive model comprising 30 cells of 5 mm length was constructed. Each cell initially contained the mineral volumes measured by X-ray diffraction (XRD) and porosity measured by (V)SANS in unreacted portions of the caprock. The fluid chemistry occupying the boundary reservoir cell was based on analyses of the modern reservoir fluids 43 with a redox state defined using the SO 4 2 À /H 2 S redox couple (sample DFS00413 with SO 4 20.7 mmol and H 2 S 0.5 mmol). The initial caprock pore fluid chemistry was calculated similar to that in equilibrium with the caprock mineralogy (reactive minerals haematite, Mg-dolomite, K-feldspar and illite), with pCO 2 , pO 2 and salinity typical of Jurassic marine shales 43 (Supplementary Table 4). Models were run assuming local fluid-mineral equilibrium and pyrite, an Fe-bearing dolomite (Mg 0.9 Fe 0.1 CaCO 3 ), and illite were allowed to precipitate (Supplementary Table 2). A D e value of 5 Â 10 À 12 m 2 s À 1 was used for all aqueous species. The model timescale was 125,000 years, with a time step of 7 days required to reach convergence of the numerical solution.
The models reproduce the key observations (Fig. 4a). The haematite dissolution front migrates B7 cm in 410 4 years. Dolomite dissolves downstream of the front and Fe-dolomite is precipitated upstream of the reaction front ( Supplementary  Fig. 2B). Pyrite is precipitated initially close to the caprockreservoir contact and migrates very slowly. K-feldspar dissolves  (4) given the time for the reaction front to develop (t 0 ; equation (3)) and the displacement distance of the front, l. Calculations are performed on the basis of Fe in equation (1). Shaded area denotes solutions for 37oqo228 m À 1 , which gives a haematite dissolution rate, k R , between 1 Â 10 À 10 and 5 Â 10 À 9 mol m À 2 s À 1 from equation (5) for a diffusion coefficient of 5 Â 10 À 12 m 2 s À 1 . Solid line shows median solution with D e ¼ 5 Â 10 À 12 m 2 s À 1 and DC 0 ¼ 2.9 mol m À 3 . Dashed lines show extreme solutions with the minimum estimate of diffusivity of 10 À 12 m 2 s À 1 and minimum DC 0 of 1.9 mol m À 3 giving maximum reaction times and a high estimate of diffusivity of 1 Â 10 À 11 m 2 s À 1 , combined with a maximum estimate of DC 0 of 3.95 mol m À 3 giving minimum reaction times. Uncertainties in the duration of caprock alteration propagated from the other variables are small compared with those from the potential range of caprock diffusivities and reaction stoichiometry.
only at the inlet, driving local precipitation of illite. The additional complexities observed with multiple zones of pyrite precipitation probably result from the episodic release of CO 2 -charged brines by the fault system 26 . Models run with two episodes of CO 2 injection reproduce a similar pattern in which the diffusion of mineral profiles into the caprock preserves a record of the evolution of reservoir fluid compositions (Fig. 4b) apparent as double peaks in pyrite and dolomite.

Discussion
This study provides evidence for the slow penetration rate of CO 2 -promoted fluid-rock reactions into clay-rich caprocks and the development of a stable reaction front without a significant positive feedback loop between porosity enhancement and the velocity of the front. SANS methods are effective at resolving micrometre-scale changes in pore network characteristics and diffusive transport properties, with the magnitude of the change of the latter being approximately half an order of magnitude. This provides evidence to refute observations from laboratory experiments, which posit a significant increase in the transport properties of caprocks following exposure to CO 2 -charged brines 14 , and CO 2 -charged brines and H 2 S 17 . Multicomponent numerical modelling of the reactive-diffusive transport reproduced the major mineralogical changes accompanying diffusion of the CO 2 , although detailed variations reflect the complex heterogeneity of natural rocks, and the cyclic variation in CO 2 charging of the underlying reservoir. The important conclusion from the analytical and numerical modelling is that the timescale for the development of the 7-cm reaction front is B10,000 to 100,000 years, consistent with the geological evidence for the duration of the CO 2 supply to the Little Grand Fault system 25 . Previous numerical models of reactive-diffusive transport of CO 2 in clay-rich caprocks, unconstrained by observations from natural systems, have predicted the development of mineral alteration profiles over length scales of metres on a 10,000 year timescale 5,9,10,12 . In addition, these modelling studies did not consider redox reactions, which as demonstrated here, are an important sink for CO 2 invading the caprock and contribute significantly to retardation of the CO 2 -diffusive transport. The results demonstrate that numerical models of reactive transport in low-diffusivity media, where mineral reaction kinetics do not limit reaction rates, successfully predict reactions and reaction rates, providing the full chemical complexity is modelled. The results also show that the mineralogy of caprocks and chemistry of potential CO 2 -charged brines will need to be considered on a case-by-case basis. Fluid-rock reaction has substantially retarded diffusive transport of the CO 2 , which would have penetrated B8 m in 100,000 years in the absence of reaction, based on a diffusion coefficient of 5 Â 10 À 12 m 2 s À 1 . Rock-buffering reactions are effective at establishing fluid-mineral equilibrium over centimetre length scales. In this caprock, CO 2 -charged fluid-mineral reactions act to enhance caprock integrity over a time period comparable to that needed for effective geological carbon storage.

Methods
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 B3-5 mm resolution, using a steel blade rock saw. Approximately 2-5 g of air-dried sample was powdered to r100 mm 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 fault 44 , located B35 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). 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°2y recorded from 2°to 92°2y. 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 models 46 . 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.

XRD measurements. XRD measurements to determine mineralogy (Supplementary
Stable isotopes. Carbonate mineral d 13 C and d 18 O were determined at the University of Cambridge on duplicate subsamples of 300-500 mg 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, (2s). d 13 C and d 18 O results are reported, relative to the VPDP and VSMOW standards, as averages of these duplicate analyses, with error bars that are the 2s s.e. of duplicate averages (Supplementary Table 2 and Supplementary Fig. 1).  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 fluidmineral equilibrium and a constant D e value of 5 Â 10 À 12 m 2 s À 1 was used for all aqueous species. The model timescale was 125,000 years with a time step of 7 days. (b) One-dimensional model with the same starting conditions as a but with two 25,000 year phases of CO 2 -saturated brine in the basal cell with an intervening 75,000 year period where the basal cell contained a CO 2 -poor brine. The results of the two episodes of CO 2 saturation are seen as double peaks in pyrite and dolomite modes as the CO 2 from each pulse diffuses into and reacts with the caprock.
Sr-isotopes. 87 Sr/ 86 Sr 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 (1s) on 158 separate measurements made during the course of these analyses. Strontium blanks were always o250 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 mm) 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 staticvolumetric 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 mm 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 o0.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 nm 2 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 B40 Å. Repeated measurements on standard materials demonstrated that the accuracy and precision of the methods described above is better than 2% (2s).
Small angle neutron scattering. Sample porosity and specific surface area was analysed using SANS and VSANS techniques (Fig. 3, Supplementary Figs 4 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 mm 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.

and 6, and Supplementary
For (V)SANS measurements, a collimated neutron beam is elastically scattered by the sample 48,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 y by Q ¼ (4p/l)sin(y/2), where l is the wavelength of the neutrons. Thus, the size range of features accessible with neutron scattering depends on the neutron wavelength l and the range in the scattering angle y. SANS data at KWS-1 were collected at wavelengths of l ¼ 0.69 nm with a wavelength distribution of the velocity selector Dl/l ¼ 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 6 Li glass scintillation detector with an active area of 60 Â 60 cm 2 . VSANS data at KWS-3 were collected at l ¼ 1.28 nm, Dl/l ¼ 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 6 Li 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 mm (rE2.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 f is obtained directly 48,50 : We assumed that the scattering intensity of pore features is directly proportional to the scattering contrast between matrix and pore 48 : Here, p Ã 1 and p Ã 2 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 fractals 48 . The SLD Dp* of each mineral was calculated as: where b i and M i are scattering length and atomic mass of the i th element in the mineral, d i is the mass density of the i th mineral and N A is the Avogadro number. The terms f 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 as 48,51 : where r 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 PRINSAS 51 . 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 (m 2 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, f, by where t 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 three 52 , 2rD v r3. 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 rQr10 À 2 Å À 1 .
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms12268 ARTICLE 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 (m 3 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 (m 2 ), DL s is the sample length (m) and DP is the pressure drop along the sample (Pa). The equivalent permeability term (k) is given by where m 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 m 2 with an average value of 9.1 Â 10 À 22 m 2 (1s ¼ 2.3 Â 10 À 22 m 2 ).
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 Stumm 41 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, (o 0 f) where o 0 is fluid velocity and f porosity, through the caprock are estimated to be o6 Â 10 À 12 m s À 1 given permeabilities of B10 À 21 m 2 and a piezometric pressure gradient of B10 4 to 4 Â 10 6 Pa m À 1 taking the hydraulic head in the Navajo 54 either across the whole Carmel Formation or the 16 cm basal seal. With effective diffusivities (D e ) of components in the fluid phase of B10 À 11 to 10 À 12 m 2 , the Péclet number (o 0 fh/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 (m 2 s À 1 ), porosity f, tortuosity, t 2 , mineral reaction rate k f (m s À 1 ) and surface area a (m 2 m À 3 ) can be written as 39 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, f 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 fD/t 2 and C 0 is solute concentration in the input fluid. Substituting equation (21) into equations (19) and (20) gives @C 0 @t 0 ¼ @ 2 C 0 @x 02 þ N D C 0 ð22Þ @f 0 S @t 0 ¼ À N D C 0 ð23Þ 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 m 2 s À 1 , Damköhler numbers are in the range 5 Â 10 7 to 9 Â 10 9 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 Â 10 3 to 3 Â 10 7 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 t 0 (s) given by or in dimensionless time, At times greater than t 0 , the reactive mineral is exhausted over distances xrl.
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 valid 39 . Solutions for fluid compositions and solid mineral modes, assuming mineral surface area remains constant, following ref. 39, in dimensionless coordinates for t 0 40, are (note l 0 ¼ 1) durations of each pulse were chosen to mimic the timescale of periodic degassing of CO 2 previously document for the local fault systems 26 .
Data availability. The authors declare that the data supporting the findings of this study are available within the article and its Supplementary Information files.