This is an open-access article distributed under the terms of the Creative Commons Attribution CC-BY 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.Open access: Article available to all readers online. This article is CCBY.


For the purpose of geological carbon storage, it is necessary to understand the long-term effects of introducing CO2 and sulfur-species into saline aquifers. CO2 stripped from the flue gas during the carbon capture process may contain trace SO2 and H2S and it may be economically beneficial to inject S-bearing CO2 rather than costly purified CO2. Furthermore, reactions between the S-CO2-bearing formation brines and formation minerals will increase pH and promote further dissolution and precipitation reactions. To investigate this we model reactions in a natural analog where CO2- and SO4-H2S bearing fluids have reacted with clay-rich siltstones. In the Mid-Jurassic Carmel formation in a cap rock to a natural CO2-bearing reservoir at Green River, Utah, a 3.1 mm wide bleached alteration zone is observed at the uppermost contact between a primary gypsum bed and red siltstone. Gypsum at the contact is ~1 mm thick and shows elongate fibers perpendicular to the siltstone surface, suggesting fluid flow along the contact. Mineralogical concentrations, analyzed by Quantitative Evaluation of Minerals by SCANning electron microscopy (QEMSCAN), show the altered siltstone region comprises two main zones: a 0.8 mm wide, hematite-poor, dolomite-poor, and illite-rich region adjacent to the gypsum bed; and a 2.3 mm wide, hematite-poor, dolomite-poor, and illite-poor region adjacent to the hematite alteration front. A one-component analytical solution to reactive-diffusive transport for the bleached zone implies it took less than 20 yr to form before the fluid self-sealed, and that literature hematite dissolution rates between 10–8 and 10–7 mol/m2/s are valid for likely diffusivities. Multi-component reactive-diffusive transport equilibrium modeling for the full phase assemblage, conducted with PHREEQC, suggests dissolution of hematite and dolomite and precipitation of illite over similar short timescales. Reaction progress with CO2-bearing, SO4-rich, and minor H2S-bearing fluids is shown to be much faster than with CO2-poor, SO4-rich with minor H2S-bearing fluids. The substantial buffering capacity of mineral reactions demonstrated by the S- and CO2-related alteration of hematite-bearing siltstones at the Green River CO2 accumulation implies that corrosion of such a cap rock are, at worst, comparable to the 10 000 yr timescales needed for carbon storage.


Storing anthropogenic carbon dioxide (CO2) in geological reservoirs is a primary option to mitigate climate change. One of the perceived risks of carbon storage is the leakage of stored carbon dioxide (CO2) through cap rocks, faults, and fractures. To quantify these risks through the life of an individual carbon storage site it is essential to understand the geochemical behavior and evolution of CO2 in geological reservoirs over a range of timescales. Best practice guidelines suggest leakage rates of less than 0.01% of the volume stored per year and storage times on the order of several thousand to 10 000 yr are required to ensure that storage sites meet the ultimate aim of mitigating climate change (Chadwick et al. 2008).

CO2 stripped from the flue gas during the carbon capture process is likely to contain trace (<5%) co-containments, i.e., SO2, H2S, NOx, and O (Talman 2015). To purify the gas further is a costly process and hence, if CO2 can be safely stored with these co-contaminants it will be economically beneficial. When SO2 dissolves in groundwater it may form H2SO3, H2SO4, or H2S depending upon redox conditions, and these compounds are stronger acids than dissolved carbonate species formed through CO2 dissolution. The resulting decreases in pH and pe are likely to result in enhanced mineral dissolution and precipitation at the gas-groundwater interface (e.g., Knauss et al. 2005; Xu et al. 2007; Wilke et al. 2012). Only a few studies have investigated the short-term impact of these coupled reactions through either batch reactor experiments (Dawson et al. 2015; Pearce et al. 2015a, 2015b) or field experiments (Bachu et al. 2005; Kaszuba et al. 2011). Over geological time even less is known about the impact of the dissolved sulfur species and CO2 on cap rocks and reservoir rocks. These long exposure times may lead to porositypermeability changes from mineral precipitation and dissolution. The study of natural analogs with S- and CO2-bearing reservoir fluids are needed to constrain long-term predictions of reservoir and cap rock alteration during CO2 storage.

Natural CO2 and CO2-SO4-bearing fluids have been migrating up two fault zones for >400 000 yr at Green River, Utah, through sandstones, mudstones, and siltstones (Burnside et al. 2013). These lithologies are typical of those proposed as geological repositories for anthropogenic CO2. Scientific drilling adjacent to the Little Grand Wash fault in 2012 has allowed study of CO2 transport mechanisms and long-term reactions between S-bearing, CO2-saturated brines and formation minerals in shallow aquifers (Fig. 1; Kampman et al. 2013, 2014, 2016; Busch et al. 2014). Drill-core collected in this 2012 drilling campaign exhibits bleaching in the reservoir sandstone, basal sections of the cap rocks related to these migrating fluids. The basal ~10 cm of a siltstone bed in the middle of the Entrada Sandstone and a claystone bed at the base of the Carmel Formation are bleached where in contact with CO2- and S-bearing reservoir fluids (see Busch et al. 2014 and Kampman et al. 2016 for further details). Small-scale bleaching is also observed within the Carmel Formation adjacent to horizontal and sub-horizontal gypsum-filled fractures. Figure 2 shows the contact between the top of a remobilized gypsum bed and the base of a red siltstone bed. This remobilized gypsum bed associated with the bleaching has a 1 mm wide fibrous overgrowth perpendicular to bedding where in contact with the siltstone, suggesting fluid flow along a previously open fracture. Numerous gypsum-filled fractures have been comprehensively documented and investigated in the core indicating gypsum-saturated fluids were common (Chen et al. 2016). Additionally, celestine and gypsum have been found as fracture fill at surface in the study area (Dockrill 2006). This is supported by the modern day fluids, which are mildly reducing, and CO2-saturated at depth (~320 m) with minor S and CH4 (Kampman et al. 2014).

In this paper we use mapping of quantitative elemental mineralogy by QEMSCAN techniques to calibrate reactive transport modeling. The modeling is used to determine the duration of the bleaching event and the sensitivity of the bleaching reactions to fluid composition; specifically, the impact of CO2-bearing fluids on clay-rich siltstones.

Regional geology

Near the town of Green River, natural CO2-bearing fluids are hosted in a series of fault bounded Jurassic sandstone reservoirs in the footwall blocks of the southerly dipping Little Grand Wash and Salt Wash normal faults, where they intersect the apex of the Green River anticline (e.g., Heath et al. 2009; Dockrill and Shipton 2010; Kampman et al. 2009, 2012, 2013; Shipton et al. 2004; Wigley et al. 2012). The CO2-bearing fluids escape to the surface through the fault damage zones of the faults as well as several abandoned petroleum exploration and water wells (Fig. 1) where they form a series of CO2 springs and geysers (Heath 2004; Kampman et al. 2009; Dockrill and Shipton 2010; Shipton et al. 2004, 2005). Ancient travertine deposits along the faults attest to CO2 leakage for at least 400 000 yr (Burnside et al. 2013). The intermittent travertine deposition is driven by the CO2-leakage with an episodicity controlled by climate driven changes in the hydraulic behavior of the faults (Kampman et al. 2012). Meteoric fluid flows in the sandstone reservoirs from recharge regions in the northern San Rafael Swell to discharge in the Green River, south of the Green River anticline (Hood and Patterson 1984). Fluid flow is parallel to the faults (west to east) where they are sealing, and toward the southeast where they are transmissive toward the fault tips. Artesian conditions prevail throughout the transmissive formations of the Palaeozoic stratigraphy, driving fluid migration up the normal faults to discharge at CO2-springs at the surface. Details of the local hydrology are discussed in Kampman et al. (2009) and the regional hydrology is discussed in Hood and Patterson (1984).

The regional burial diagenesis of the Jurassic strata includes development of hematite and goethite grain coatings, giving the rock its characteristic red color (Cullers 1995; Trimble and Doelling 1978). Proximal to hydrocarbon seeps and CO2 accumulations, extensive bleaching has occurred (e.g., Beitler et al. 2003, 2005; Parry et al. 2004; Wigley et al. 2012). At Green River, Wigley et al. (2012) document bleaching at the base of the stratigraphy in Salt Wash Graben, ascertaining that the reducing fluid (i.e., CO2-S-bearing brines) must be more dense than the formation fluid. The modern fluids in the Entrada and Navajo Sandstone are CO2- and SO4-rich with a low pH of 5.1. Furthermore, Eh measurements of spring and geyser waters indicate the fluids are mildly reducing (0 to –50 mV; Kampman et al. 2014). Analyses of fluid inclusions by Wigley et al. (2012) and exsolved gas compositions from the springs by Kampman et al. (2014) show these fluids also contain trace quantities of CH4 (0–28%), a known reducing agent. Kampman et al. (2014, 2016) inferred that reduced sulfur must be playing a role in the bleaching given the mildly reducing groundwater (0 to –50 mV that implies that the fluids with ~20 mmol/L SO4 contain ~0.5 mmol/L H2S) coupled with the evidence for precipitation of pyrite in the basal portions of the Entrada and Carmel cap rocks. CO2-bearing core experiments conducted by Purser et al. (2014) on red Entrada Sandstone showed H2S is most likely the main driver for bleaching. Using 1 mg/L of thioacetamide (which breaks down on the addition to water releasing reduced sulfur) the authors were able to reduce all the Fe3+ in the core sample to Fe2+ (mobilizing 15 mg/L of iron), thus bleaching the entire sample in the 6 month experiment. Whereas, using 5% CH4 they were only able to mobilize at best one fifth (3 mg/L) of the iron in 6 months. It is evident that the kinetics of the hematite dissolution reaction are much faster with H2S rather than with CH4 as the reducing agent. Given this evidence, H2S is considered to be the primary reductant in the Green River CO2 accumulation.

Fluid migration and subsequent alteration in the region likely occurred either: (1) pre-CO2 accumulation during uplift associated with the formation of the Colorado Plateau after the Eocene (S-bearing, CO2-poor fluids), or (2) synchronous with the CO2 accumulation and possibly relating to the periodic recharging of the reservoir during the Quaternary (S-bearing, CO2-saturated fluids; Kampman et al. 2012).

Local geology

Scientific drilling into the Green River CO2 accumulation (see Kampman et al. 2013) intersected the Entrada Sandstone, Carmel Formation and the Navajo Sandstone. Within the Carmel Formation a red siltstone layer is bleached where it is in contact with a gypsum vein adjacent to the top of a primary gypsum bed. The Carmel Formation is a 50 m thick complex lithologically heterogeneous package representing units that were deposited in environments ranging from aeolian, fluvial, coastal sabkha to shallow and open marine (Blakey et al. 1983, 1988). The documented bleached contact occurs at 158.34 m below surface within a package of interbedded, unfossiliferous red and gray siltstone/shale and bedded gypsum/anhydride (Fig. 2). Extensional faulting within the Carmel Formation has resulted in the re-mobilization of many of these gypsum and anhydride beds, along with the brecciation of siltstone beds and formation of gypsum-filled open fractures. It is above one of these gypsum filled open fractures adjacent to the top of a primary gypsum bed where 3.1 mm of a red siltstone bed is bleached. The bleaching is suggestive of penetration of reducing CO2-saturated fluids by diffusion and subsequent alteration of the siltstone mineralogy.

Analytical methods

The small scale of this alteration front means using standard techniques (i.e., subsampling for XRD and XRF measurements) for determining quantitative miner-alogical and chemical changes is challenging. As such, in this study we take a different approach by quantifying mineralogical profiles from QEMSCAN analyses and coupling these with reactive transport modeling. QEMSCAN is a proven technique recently becoming popular with industry and research organizations including studies relating to CO2 storage (e.g., Petrel Sub-Basin, Northern Territory, Australia; Consoli et al. 2013). Generally good agreement has been seen with XRD analyses (Wunsch et al. 2014; Farquhar et al. 2015; Pearce et al. 2016).

Mineralogical and elemental maps were made of a thin section across a gypsumsiltstone reaction front located in the drill-core at 158.34 m below surface from the drill-hole CO2W55 (Fig. 2; Kampman et al. 2013). The quantitative elemental maps were obtained using a Quanta 650F, field emission gun (FEG) scanning electron microscope (SEM), equipped with two Bruker XFlash 6130 energy dispersive spectrometers (EDS) at the Department of Earth Sciences, University of Cambridge. The fully automated system includes an automated spectrum acquisition and classification procedure. Analyses were performed by obtaining field-scans (i.e., phase maps) across the areas of interest thus providing a complete characterization of particle surfaces above a predefined electron backscatter threshold. The brightness coefficients were calibrated against quartz, gold, and copper. Spectra were collected at 15 kV and 10 nA with 2000 total X-ray counts per pixel at a 1.5 μm spacing, and compared to a Species Identification Protocol (SIP) that discriminates minerals on the basis of their characteristic X-ray and electron backscatter intensities. Data was processed further through a FEI software package under development, NanoMin, which has a better capacity in dealing with solid solutions allowing for superior resolving power of the chemistry in each pixel. This software package was used to obtain mineralogical and elemental maps illustrating the textural and chemical relationships across the reaction front. XRD is better for clay identification, however it requires a larger amount of sample and only gives a bulk measurement. In comparison, QEMSCAN is particularly powerful for small samples. Additionally, it provides information on the relationships between morphology, texture, and chemistry in a mapped area, which is useful for this application.

Mineralogical changes

The early diagenetic mineralogy of the siltstones in the Carmel Formation is dominated by illite, quartz, dolomite, ankerite, K-feldspar, albite, and hematite (Appendix1 Table 1; Fig. 3). The contacts between the gypsum/anhydride beds and the siltstone beds mapped in the drillcore (Kampman et al. 2014) show no evidence of alteration with the exception of one upper contact of a gypsum bed at 158.34 m where a ~3.1 mm thick zone of alteration occurs in the siltstone adjacent to the gypsum bed. Chen et al. (2016) measured the hydration waters of this gypsum bed at 158.48 m, the sample comprised 19.7% H2O. Figures 3 and 4 show the major mineralogical profiles of the alteration zone adjacent to the gypsum bed as revealed by QEMSCAN analyses. The mineralogy across the reaction front is determined by the NanoMin software using mineral compositions of:

  • Albite: NaAlSi3O8

  • Ankerite: CaFe0.7Mg0.3(CO3)2

  • Dolomite: CaMg(CO3)2

  • Hematite: Fe2O3

  • Illite: K0.97(H3O)0.4Al2.46Fe0.3Mg0.17Ti0.03Si3.07O7.15(OH)2·(H2O)

  • Orthoclase/Microcline: K(AlSi3O8)

  • Quartz: SiO2

It is evident from Figures 3 and 4 that there are two dominant alteration zones: (1) a ~2 mm wide zone of Fe dissolution adjacent to the sharp bleached-red transition, and (2) a ~1 mm zone of Fe, dolomite, orthoclase/microcline, and albite dissolution, and the precipitation of illite in the siltstone adjacent to the gypsum bed. The QEMSCAN analysis indicates that only minor hematite is present in the unaltered sample and instead the Fe is bound in Fe-illite and ankerite (ankerite has a similar spatial distribution to dolomite). Whether the hematite is too fine-grained to resolve or Fe-illite is actually present is not known. XRD analysis of the compositionally similar Carmel cap rocks suggests the Fe is bound in hematite (Kampman et al. 2016). For the purpose of modeling the geochemical profile it is assumed the Fe is in hematite.

An ~1 mm zone of elongate, fibrous gypsum is present at the contact between the gypsum and the siltstone with the fibers growing with the growth generally perpendicular to the contact (Fig. 5). Fibrous gypsum normally precipitates within fractures to form veins (e.g., Machel 1985; Raman and Ramdas 1954; Warren 2006), therefore the presence of these structures at the contact implies that the vein has undergone extension, and thus, at one point in time, a fracture opened between these two sedimentary units. This extensional period causing the open fracture could either be associated with the uplift of the Colorado Plateau in the Cenozoic (Flowers 2010) or with later de-glaciation events during crustal flexing from the glacial unloading of Lake Bonneville, Cutler Dam, or Little Valley (Kampman et al. 2012). The observed petrological changes in the adjacent bleached zone are consistent with a series of dissolution reactions involving CO2, SO2, and H2S that result in the removal of hematite, dolomite, orthoclase, and albite via stoichiometries such as:

Hematite dissolution:  
Dolomite dissolution:  
Orthoclase/Microcline dissolution:  
Albite dissolution:  

The released cations of Fe, Mg, K, and Al can then react form illite, i.e., by the illitization of kaolinite (i.e., Kohler et al. 2009).

An overall loss of Mg, Fe, Na, and Ca at the siltstone-gypsum bed boundary suggests the reacting fluids transported these cations out of the system. The documented petrological changes are consistent with the stoichiometry of the reaction:  

Reactive transport modeling

The redistribution of chemical components in groundwater is controlled by four principal processes: advection, mechanical dispersion, molecular diffusion, and the rate of chemical reaction (e.g., Bear 1972; Bethke 2008). Siltstones in the Carmel formation have low permeabilities, < 10–18 m2, effective solute diffusivities, De, in the fluid phase on the order of 10–11 to 10–12 m2 (Kampman et al. 2016), and a piezometric pressure gradient between 104 and 106 Pa/m (taken from the hydraulic head data of the Navajo Sandstone; Hood and Patterson 1984; Kampman et al. 2016). These parameters imply a flow rate (ω0ϕ) of <6 × 10–10 m/s for the siltstones, and Peclet numbers (Pe = ω0φh/De), the ratio of diffusive to advective timescales, between 10–7 and 10–1. As Pe << 1 diffusive transport dominates and advection can be neglected in the transport equations.

To further understand the governing processes of a reactive system we use a simplified one component analytical solution, coupling chemical reactions and physical transport processes.

One component transport model

For one-dimensional transport, perpendicular to the vein margin, the diffusion-reaction equation describing the change in solute concentration C (mol/m3) with time t (s) and distance, x (m) can be written as:  
where De is the effective diffusion coefficient, kf is the mineral reaction rate (m/s), α is the mineral surface area (m2/m3), and Ceq is the solute concentration at equilibrium. The variation in the mineral molar concentration in the rock, φS (mol/m3) with time t, as controlled by surface dissolution, is modeled by a linear kinetics as:  
Following Kampman et al. (2016), the following transformations to dimensionless variables are made:  
where the prime (′) indicates a dimensionless variable and l is an appropriate length scale. The dimensionless equations describing the changing solute concentration for one-dimensional diffusive transport with mineral reaction of over time, and the variation of mineral moles in dimensionless time (t′) are, respectively:  
The solutions for fluid compositions and hematite mineral modes are controlled by the dimensionless Damköhler number, ND, the ratio between the rate of reaction and the rate of transport of the fluid. For diffusion dominated reactive transport (Eq. 1) the Damköhler number is defined by:  

For a mildly reducing, S-bearing, CO2-saturated fluid similar to that of the downhole Navajo fluids (Kampman et al. 2014) hematite dissolution rates (kR) of 10–8 to 10–7 mol/m2/s are expected (Dos Santos Afonso and Strumm 1992). Hematite appears to be finely disseminated in the siltstone, thus grain surface areas (α) are expected to be similar to the Carmel cap rock and are estimated to be ~106 m2/m3. The transport distance, l (m), is taken as the displacement of the reaction front, 3.1 mm. Damköhler numbers for the gypsum-siltstone alteration zone in the Carmel Formation are therefore estimated to be between 2 and 15. This implies that the timescale of the system is diffusion limited and it can be assumed that the reacting fluids will be close to equilibrium with the minerals. Taking these assumptions, we use the following analytical solutions to the mass transport equations, as presented in Lichtner (1988).

After onset of hematite dissolution, the time (τ0) it takes for hematite at the vein contact to become exhausted is described, assuming a constant hematite surface area, by:  
where φS is the dimensionless volume of hematite initially present. For times, t > τ0, hematite is exhausted over distances xl.
Once this hematite dissolution front has formed (t > τ0), the time it takes for the front to migrate an observed distance, l, is calculated by:  
where q (m–1) is the exponential constant giving the length scale over which fluid returns to equilibrium with the initial mineralogy and defined as:  
Note that the Damköhler number, ND = (ql)2.
The volume fraction of hematite across the alteration front, i.e., the shape of the hematite concentration profile, over time, t, and at a distance, x, downstream of the reaction front is given by:  
Figure 6 shows the calculated position of the hematite reaction front at time, t, the geometry of which is described by Equation 7. The profile primarily reflects hematite dissolution and diagenetic variations in primary hematite concentrations across the reaction front. The geometry of this profile, assuming a constant De, is solely dependent on the kinetics of hematite dissolution. Taking an average of the hematite volume in the red unaltered region, a line of best fit across the profile is determined using the least-squares method. Exponential constant, q, in Equation 7 is therefore calculated to be 896 ± 150 m–1 (1σ), however, q could be as small as 472 m–1. With a measured reaction front distance of l = 0.0031 m, an estimated effective diffusivity, De = 5 × 10–12 m2/s, an estimated hematite surface area of 106 m2/m3 and an infiltrating fluid compositions akin to the downhole Navajo Sandstone fluid samples we can infer the Damköhler number, timescales, and dissolution rates (Fig. 7). The calculated values of q imply a Damköhler number between 2.1 and 7.7, timescales of reaction between 4 and 20 yr, for hematite dissolution rates between 2 × 10–8 and 8 × 10–6 mol/m2/s (gray shaded area; Fig. 7).

Numerical reactive transport model

Modeling the effects of multiple minerals, components, and simultaneous reactions requires the mass transport equations to be evaluated numerically. Such modeling of multiple mineral profiles (Figs. 34) provides multiple constraints. As calculated above, low Peclet numbers and high Damköhler numbers indicate the fluid transport is diffusion dominant and the fluids are close to equilibrium with the host minerals. As such, we assume instantaneous equilibrium and diffusion-only transport in the models below (to make modeling simpler and thus time efficient). The modeling of the alteration was conducted using PHREEQC 3.3.9 and the llnl.dat database (Parkhurst and Appelo 1999).

Key parameters and variables. The 3.1 mm wide gypsumsiltstone alteration zone was modeled as a 1 cm long 1D reactive-diffusive model comprising of 20 cells of 0.5 mm length. The choice of input parameters to the reactive transport model can have a significant impact on the conclusions reached. Variables, including the thermodynamic parameters of the relevant mineral and aqueous species, and their diffusivities are assumed to be constant. These variables are taken from the thermodynamic database llnl.dat, with diffusivities taken from direct measurements of the similar Carmel cap rock. The porosities and mineral volumes are based on quantitative elemental mineralogy (QEMSCAN) maps (Appendix1 Table 1). Phases allowed to dissolve and precipitate in this model were hematite, dolomite, orthoclase/microcline, albite, ankerite, gypsum, illite, and quartz. The phases were chosen based on the mineralogical maps across the reaction front. The thermodynamic parameters for ankerite [Mg0.3Fe0.7Ca(CO3)2] were taken from TOUGHREACT's thermodynamic database (Xu et al. 2011). The initial pore fluid composition of the siltstone was taken to be a fluid calculated to be in equilibrium with the unaltered siltstone mineralogy, and with pCO2, pO2, and salinity estimates as that for typical Jurassic marine shales (Turrero et al. 2006; Kampman et al. 2016; Table 1). The infiltrating fluid was taken to be a fluid in equilibrium with gypsum and saturated with CO2.

Two sets of simulations were run to investigate whether the siltstone alteration zone adjacent to the gypsum vein was associated with the infiltration of: (1) a CO2-poor, SO4-rich fluid likely relating the uplift of the Colorado Plateau after the Eocene, or (2) a CO2-SO4-saturated fluid relating to the recharge phase of the CO2-accumulation due to fluid movements during the Quarternary, perhaps triggered by deglaciations (Kampman et al. 2012). Since the composition of the bounding fluid in the vein is unknown, the basal Navajo Sandstone downhole fluid sample (DFS004; Kampman et al. 2014) is used varying its CO2 content from saturated (Model 1) to CO2-poor (Model 2). The starting compositions of the bounding fluids are shown in Table 1. A third set of models (Model 3) were run to test the sensitivity of the fronts to effective diffusivity, with the CO2-SO4-saturated fluid as the infiltrating fluid. All infiltrating fluids were initially equilibrated with gypsum. The models used PHREEQC assuming diffusive transport in the fluid phase and local mineral-fluid equilibrium and were run for 100 yr with a time step of 7 days.

Uncertainties. A Monte Carlo approach is used to assess uncertainties contributions from infiltrating fluid compositions and effective diffusivity. Additional uncertainties in the numerical modeling are associated with thermodynamic databases and the assumption of instantaneous fluid-rock equilibrium.

Model 1 and Model 2 were run for 1000 Monte Carlo iterations with the components Al, Ca, Mg, K, Na, Si, SO42–, H2S, Pe, and dissolved inorganic carbon (DIC; Table 1) in the bounding fluid being randomly varied in each iteration by ±10% (normal random distribution) relative to their listed concentrations. Cl concentrations were adjusted to maintain charge balance as Cl concentration are independent of reactions. In Model 3, 10 iterations were run where effective diffusivity was varied and bounding fluid composition kept constant for the CO2-saturated case only. Of the 1000 simulations in Models 1 and 2, 13 of the CO2-saturated simulations and 3 of the 1000 CO2-poor simulations fail to converge within the 100 iterations specified by default in PHREEQC's numerical equation solver. None of the 10 effective diffusivity simulations in Model 3 failed to converge. Changes in DIC, Pe, H2S, and SO42– have the largest impact on final mineral modes.


The results of the reactive transport models are shown after 6 yr of diffusion in Figures 8 and 9. A 6 yr time step is shown because it is the best fit to observations, this is discussed further below. The hematite and dolomite reaction fronts progress quicker with the CO2-SO42–-saturated infiltrating fluid. Critically the CO2-poor infiltrating fluid did not precipitate illite in any of the model runs. The peak in illite concentration in the core sample fracture is observed by quantitative elemental mineralogy (Fig. 4) adjacent to the gypsum bed, and is reproduced best by the CO2-SO42–-saturated model, suggesting the alteration is associated with the CO2-accumulation. Additionally, the velocity of the hematite reaction front is directly related to the concentration of H2S and the H2S/SO42– ratio in the infiltrating fluid. Faster front velocities are associated with a greater H2S concentration in the infiltrating fluid. In several simulations in both Model 1 and Model 2 dolomite precipitates at x = 0, where the bounding fluid is dolomite saturated. However, there is no mineralogical evidence for dolomite precipitation at the margin of the vein or in the basal contacts of siltstones with the Navajo sandstone in the Green River CO2 accumulation, which may be a consequence of kinetic limitations on dolomite precipitation rates. In the absence of diffusivity measurements, effective diffusivity for the siltstone is approximated by the stratigraphically lower Carmel claystone cap rock (5 × 10–12 m2; Kampman et al. 2016). In general, siltstones show slightly larger diffusivity than claystone. As shown in Figure 9, increasing effective diffusivity to 1 × 10–11 m2 increases the velocity of the reaction front with the doubling of diffusivity causing the expected

increase for diffusive transport. The uncertainty in diffusivity is a fundamental limitation on the inference of timescales.

Both the analytical hematite model and the multi-component PHREEQC reactive-diffusion model using the gypsum-saturated basal downhole fluid sample, give similar reproductions of the observed petrological trends. The solutions of the analytical model imply that the hematite reaction front took between 4 and 20 yr to propagate 3.1 mm given hematite reaction rates of the order of 2 × 10–8 and 8 × 10–6 mol/m2/s. A comparison at a 6 yr time step between the CO2-SO4-saturated model numerical model and petrological observations from QEMSCAN is shown in Figure 10. The model predicts dissolution of albite, hematite, and dolomite upstream of the hematite reaction front, and growth of secondary illite adjacent to the gypsum bed. The latter is not predicted by models of CO2-poor fluids. Discrepancies between the observed mineral profiles and the CO2-saturated simulations can be explained by the simplifying assumptions inherent in the models, that is, homogenous effective diffusivity, homogenous initial mineral modes and local fluid-solid equilibrium between fluids and mineral phases. A least-squares fit between the Monte Carlo simulations and the observed mineralogical profile gave a best fit at 6.1–1.0+2.0 yr (uncertainties reflecting uncertainties in fluid compositions) since the start of CO2-SO4-saturated brine injection (Fig. 10). This is comparable to the time estimates of the analytical model for hematite dissolution. The estimate of diffusivity is the main source of uncertainty where timescales as the inverse of diffusivity. If the fracture was actively precipitating gypsum during this time it would give an average linear extension rate of gypsum ~0.17 mm/yr. It is very unlikely precipitation occurred at a uniform rate though time, but this rate is comparable to the rapid carbonate extension rates reported in Frery et al. (2017) for veins in travertine mounds near Green River.


This study has shown that QEMSCAN is an extremely useful technique for investigating and quantifying the mineralogy of small scale alteration fronts. These observed mineralogical profiles (as measured by QEMSCAN) in the siltstone adjacent to the gypsum bed can be explained in terms of the interplay of mineral reactions and diffusive transport of a reactive infiltrating fluid. The alteration comprises two main zones a ~0.8 mm wide hematite-poor, dolomite-poor, illite-rich region adjacent to a gypsum bed and a ~2.3 mm wide hematite-poor, dolomite-poor, illite-poor region adjacent to the hematite alteration front. There is an apparent decrease in porosity adjacent to the gypsum-siltstone contact (0 mm) from 0 to 0.8 mm and increase in porosity from 0.8 to 1.7 mm compared to the unaltered red siltstone region. Overall there is no net change in porosity between the altered region and the unaltered region, however, without BET measurements it is difficult to draw confident conclusions about the effect the CO2-S-bearing fluids had on the sealing capabilities of the siltstone. The integration of quantitative elemental mineralogical maps and reactive transport modeling has given consistent results with an estimated reaction time of 4 to 20 yr from analytical modeling and ~6 yr from numerical modeling required to produce the observed mineralogical alteration profiles given the assumed diffusivity of 5 × 10–12 m2·s–1. It is also clear from the modeling that the illite reaction profile depends on the dissolution of orthoclase/microcline and dolomite, however, it could also be due to primary variations in clay content. Unlike the Carmel cap rock reaction profile discussed in Kampman et al. (2016), the Entrada cap rock (Busch et al. 2014), and the Entrada Sandstone reaction front around Salt Wash Graben (Wigley et al. 2012, 2013), there are no concentration peaks in metal sulfides or metal oxides upstream of the hematite reaction front. This is likely due to the lower concentration of Fe in the unaltered siltstone, with ~4 and 2 wt% hematite in the unaltered Carmel and Entrada cap rocks, respectively, compared to a hematite concentration of 0.4 wt% in this study (Busch et al. 2014; Kampman et al. 2016). At long timescales, assuming self-sealing does not occur, the numerical model predicts it would take 5000 yr for the bleaching front to propagate ~10 cm. This is two orders of magnitude faster than >100 000 yr predicted for the Carmel and Entrada cap rocks (Busch et al. 2014; Kampman et al. 2016). This discrepancy is due to a greater concentration of hematite in the unaltered Entrada and Carmel cap rocks. It is clear that the greater the hematite concentration, the stronger the buffering capacity of the CO2-S fluid-rock reactions will be.

Other studies of CO2 altered siltstone and claystone caprocks have found similar mineralogical trends to those seen in this study (i.e., Busch et al. 2014; Kampman et al. 2016). Gypsum precipitation has been documented in the Carmel caprock where gypsum fills bedding parallel fractures (Kampman et al. 2016). Kohler et al. (2009) found that at low temperature and pressures in clayey caprocks CO2 alteration resulted in the formation of Fe2+- and K+-enriched illites. Overall this study agrees with those of Busch et al. (2014) and Kampman et al. (2016) who found that CO2-S-related fluid-rock reactions have the capability to buffer the corrosive nature of the fluids, allowing cap rocks to maintain their sealing capabilities.


Our results are important for geological carbon storage as they show that the addition of sulfur to CO2 producing potentially corrosive fluids will not degrade the sealing capabilities of hematite-bearing mud- and siltstone-rich cap rocks in carbon storage sites over the lifetime of the storage projects. Furthermore, the timescales calculated from the modeling show that, even adjacent to an active fault zone, fluid-filled cracks are sealed within a few years, further impeding migration of CO2-charged fluids. This coupling of high spatial resolution quantitative elemental and phase mapping by QEMSCAN with thermodynamic modeling should also have a wider application to the study small-scale alteration fronts in several different geological settings, particularly in stockwork veining and associated alteration in economic mineral deposits, especially since standard methods require larger sample volumes, or may miss small-scale variability and so only provide phase information or bulk compositions.


We thank Shell Global Solutions for giving us access to the samples and FEI for allowing us access to the proprietary software NanoMin. Funding for Carbon Storage research at the University of Cambridge was provided by the Natural Environment Research Council (NERC) to the CRIUS consortium (NE/F004699/1), Shell Global Solutions International, and the U.K. Department of Energy and Climate Change (DECC) through a CCS Innovation grant. A.M. was supported by an EPSRC doctoral training grant.


Deposit item AM-18-26138 Appendix Table 1. Deposit items are free to all readers and found on the MSA web site, via the specific issue's Table of Contents (go to