A siltstone reaction front related to CO 2 - and sulfur-bearing fluids: Integrating quantitative elemental mapping with reactive transport modeling k

For the purpose of geological carbon storage, it is necessary to understand the long-term effects of introducing CO 2 and sulfur-species into saline aquifers. CO 2 stripped from the flue gas during the carbon capture process may contain trace SO 2 and H 2 S and it may be economically beneficial to inject S-bearing CO 2 rather than costly purified CO 2 . Furthermore, reactions between the S-CO 2 -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 CO 2 - and SO 4 -H 2 S bearing fluids have reacted with clay-rich siltstones. In the Mid-Jurassic Carmel formation in a cap rock to a natural CO 2 -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/m 2 /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 CO 2 -bearing, SO 4 -rich, and minor H 2 S-bearing fluids is shown to be much faster than with CO 2 -poor, SO 4 - rich with minor H 2 S-bearing fluids. The substantial buffering capacity of mineral reactions demonstrated by the S-and CO 2 -related alteration of hematite-bearing siltstones at the Green River CO 2 accumulation implies that corrosion of such a cap rock are, at worst, comparable to the 10 000 yr timescales needed for carbon storage.


IntroductIon
Storing anthropogenic carbon dioxide (CO 2 ) 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 (CO 2 ) 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 CO 2 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).
CO 2 stripped from the flue gas during the carbon capture process is likely to contain trace (<5%) co-containments, i.e., SO 2 , H 2 S, NOx, and O (Talman 2015).To purify the gas further is a costly process and hence, if CO 2 can be safely stored with these co-contaminants it will be economically beneficial.When SO 2 dissolves in groundwater it may form H 2 SO 3 , H 2 SO 4 , or H 2 S depending upon redox conditions, and these compounds are stronger acids than dissolved carbonate species formed through CO 2 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. 2015aPearce et al. , 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 CO 2 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 CO 2 -bearing reservoir fluids are needed to constrain long-term predictions of reservoir and cap rock alteration during CO 2 storage.
Natural CO 2 and CO 2 -SO 4 -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 CO 2 .Scientific drilling adjacent to the Little Grand Wash fault in 2012 has allowed study of CO 2 transport mechanisms and long-term reactions between S-bearing, CO 2 -saturated brines and formation minerals in shallow aquifers (Fig. 1; Kampman et al. 2013Kampman et al. , 2014Kampman et al. , 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 CO 2 -and S-bearing reservoir fluids (see Busch et al. 2014 andKampman 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 CO 2 -saturated at depth (~320 m) with minor S and CH 4 (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 CO 2 -bearing fluids on clay-rich siltstones.

regIonAl geology
Near the town of Green River, natural CO 2 -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. 2009Kampman et al. , 2012Kampman et al. , 2013;;Shipton et al. 2004;Wigley et al. 2012).The CO 2 -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 CO 2 springs and geysers (Heath 2004;Kampman et al. 2009;Dockrill and Shipton 2010;Shipton et al.FIgure 1. Geological map of the Green River anticline showing locations of the Little Grand Wash and Salt Wash Graben normal fault systems, CO 2 -springs and location of drill-hole CO2W55.Structure contours are for the top surface of the Navajo Sandstone, the main CO 2bearing reservoir (modified after Kampman et al. 2013).
FIgure 2. Sedimentary log, core photos, and plane-polarized light (PPL) images of the bleached alteration zone at the upper contact of the gypsum bed.Images show the bleached siltstone reaction zone at the upper contact of a gypsum bed.2004,2005).Ancient travertine deposits along the faults attest to CO 2 leakage for at least 400 000 yr (Burnside et al. 2013).The intermittent travertine deposition is driven by the CO 2 -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 CO 2 -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 CO 2 accumulations, extensive bleaching has occurred (e.g., Beitler et al. 2003Beitler et al. , 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., CO 2 -S-bearing brines) must be more dense than the formation fluid.The modern fluids in the Entrada and Navajo Sandstone are CO 2 -and SO 4 -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 CH 4 (0-28%), a known reducing agent.Kampman et al. (2014Kampman et al. ( , 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 SO 4 contain ~0.5 mmol/L H 2 S) coupled with the evidence for precipitation of pyrite in the basal portions of the Entrada and Carmel cap rocks.CO 2 -bearing core experiments conducted by Purser et al. (2014) on red Entrada Sandstone showed H 2 S 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 Fe 3+ in the core sample to Fe 2+ (mobilizing 15 mg/L of iron), thus bleaching the entire sample in the 6 month experiment.Whereas, using 5% CH 4 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 H 2 S rather than with CH 4 as the reducing agent.Given this evidence, H 2 S is considered to be the primary reductant in the Green River CO 2 accumulation.
Fluid migration and subsequent alteration in the region likely occurred either: (1) pre-CO 2 accumulation during uplift associated with the formation of the Colorado Plateau after the Eocene (S-bearing, CO 2 -poor fluids), or (2) synchronous with the CO 2 accumulation and possibly relating to the periodic recharging of the reservoir during the Quaternary (S-bearing, CO 2 -saturated fluids; Kampman et al. 2012).

locAl geology
Scientific drilling into the Green River CO 2 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(Blakey et al. , 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 CO 2 -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 mineralogical 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 CO 2 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 (Appendix 1 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)  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 CO 2 , SO 2 , and H 2 S that result in the removal of hematite, dolomite, orthoclase, and albite via stoichiometries such as: Orthoclase/Microcline 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).

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 m 2 , effective solute diffusivities, D e , in the fluid phase on the order of 10 -11 to 10 -12 m 2 (Kampman et al. 2016), and a piezometric pressure gradient between 10 4 and 10 6 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 (w 0 j) of <6 × 10 -10 m/s for the siltstones, and Peclet numbers (Pe = w 0 jh/D e ), 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/m 3 ) with time t (s) and distance, x (m) can be written as: where D e is the effective diffusion coefficient, k f is the mineral reaction rate (m/s), a is the mineral surface area (m 2 /m 3 ), and C eq is the solute concentration at equilibrium.The variation in the mineral molar concentration in the rock, j S (mol/m 3 ) 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:

∂C
The solutions for fluid compositions and hematite mineral modes are controlled by the dimensionless Damköhler number, N D , 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, CO 2 -saturated fluid similar to that of the downhole Navajo fluids (Kampman et al. 2014) hematite dissolution rates (k R ) of 10 -8 to 10 -7 mol/m 2 /s are expected (Dos Santos Afonso and Strumm 1992).Hematite appears to be finely disseminated in the siltstone, thus grain surface areas (a) are expected to be similar to the Carmel cap rock and are estimated to be ~10 6 m 2 /m 3 .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 (t 0 ) it takes for hematite at the vein contact to become exhausted is described, assuming a constant hematite surface area, by: where j S ∞ is the dimensionless volume of hematite initially present.For times, t > t 0 , hematite is exhausted over distances x ≤ l.
Once this hematite dissolution front has formed (t > 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, N D = (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 D e , 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 leastsquares method.Exponential constant, q, in Equation 7 is therefore calculated to be 896 ± 150 m -1 (1s), 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, D e = 5 × 10 -12 m 2 /s, an estimated hematite surface area of 10 6 m 2 /m 3 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/m 2 /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.3-4) 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.9and the llnl.datdatabase (Parkhurst and Appelo 1999).
Key parameters and variables.The 3.1 mm wide gypsumsiltstone alteration zone was modeled as a 1 cm long 1D reactivediffusive 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 (Appendix 1 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 [Mg 0.3 Fe 0.7 Ca(CO 3 ) 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 pCO 2 , pO 2 , 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 CO 2 .
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 CO 2 -poor, SO 4 -rich fluid likely relating the uplift of the Colorado Plateau after the Eocene, or (2) a CO 2 -SO 4 -saturated fluid relating to the recharge phase of the CO 2 -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 CO 2 content from saturated (Model 1) to CO 2 -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 CO 2 -SO 4 -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, SO 4 2-, H 2 S, 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 varied and bounding fluid composition kept constant for the CO 2 -saturated case only.Of the 1000 simulations in Models 1 and 2, 13 of the CO 2 -saturated simulations and 3 of the 1000 CO 2 -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, H 2 S, and SO 4 2-have the largest impact on final mineral modes.

results
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 CO 2 -SO 4 2--saturated infiltrating fluid.Critically the CO 2 -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 CO 2 -SO 4 2--saturated model, suggesting the alteration is associated with the CO 2 -accumulation.Additionally, the velocity of the hematite reaction front is directly related to the concentration of H 2 S and the H 2 S/SO 4 2-ratio in the infiltrating fluid.Faster front velocities are associated with a greater H 2 S 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 FIgure 6. Volume fraction of hematite across the gypsum-siltstone reaction front.Squares represent the observed lithological changes as determined by quantitative mapping.The solid curve is least-squares best fit of the hematite profile to Equation 7 with the shaded area representing 1s St.dev.The initial hematite volume, j S ∞ , is the average volume over the red siltstone region.Dashed line shows profile for the minimum limit for q values of 427 m -1 .5given the time for the reaction front to develop (t 0 ; Eq. 4).The shaded area represents the solutions for 427 < q < 896 m -1 , implying timescales of reaction between 4 and 20 yr and hematite dissolution rates between 1 × 10 -9 and 5 × 10  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 CO 2 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 m 2 ; 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 m 2 increases the velocity of the reaction front with the doubling of diffusivity causing the expected 2 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/m 2 /s.A comparison at a 6 yr time step between the CO 2 -SO 4 -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 CO 2 -poor fluids.Discrepancies between the observed mineral profiles and the CO 2 -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 CO 2 -SO 4 -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.

dIscussIon
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,  1.The black lines are the mean of all these models.The precipitation of illite in the CO 2 -saturated model matches the observed mineralogical profiles the best, suggesting the fluids responsible for the alteration were CO 2 -saturated.
FIgure 9. Diffusion-reaction modeling of gypsum-siltstone alteration zone in PHREEQC at with a CO 2 -SO 4 -saturated infiltrating fluid (see Table 1 for composition).1D models were run exactly like the models presented in Figure 8, except the effective diffusivity (D e ) was varied instead of fluid composition.Three models are plotted, with D e of 1 × 10 -12 m 2 , 5 × 10 -12 m 2 , and 1 × 10 -11 m 2 , respectively.The plot shows the effect diffusivity has on the propagation speed of the reaction fronts.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 CO 2 -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 m 2 •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(Wigley et al. , 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 CO 2 -S fluid-rock reactions will be.
Other studies of CO 2 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 CO 2 alteration resulted in the formation of Fe 2+ -and K + -enriched illites.Overall this study agrees with those of Busch et al. (2014) and Kampman et al. (2016) who found that CO 2 -S-related fluid-rock reactions have the capability to buffer the corrosive nature of the fluids, allowing cap rocks to maintain their sealing capabilities.

IMPlIcAtIons
Our results are important for geological carbon storage as they show that the addition of sulfur to CO 2 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 CO 2 -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.

FIgure 3 .
FIgure 3. Images of the alteration zone that occurs between an upper contact of a gypsum bed and a siltstone bed.(a) Plane-polarized light (PPL) photomicrograph of the region that was analyzed with QEMSCAN.(b) Quantitative elemental mineralogy over the bleached/unbleached contact in the siltstone.

FIgure 4 .
FIgure 4. Images collected from QEMSCAN analyses of the alteration zone, presented are the backscatter electron (BSE) map, and quantitative elemental mineralogy maps; illite, kaolinite, dolomite, albite, orthoclase/ microcline, and hematite.FIgure 5. PPL photomicrographs of fibrous gypsum adjacent to the altered siltstone.The occurrence of fibrous gypsum at the contact implies that at one point in time there was an open fracture between these two sedimentary units allowing the flow of reactive fluids.

FIgure 7 .
FIgure 7. Plot showing the solution to the analytical one-component reactive transport equations.The duration it takes for the hematite reaction front to migrate 3.1 mm calculated from Equation5given the time for the reaction front to develop (t 0 ; Eq. 4).The shaded area represents the solutions for 427 < q < 896 m -1 , implying timescales of reaction between 4 and 20 yr and hematite dissolution rates between 1 × 10 -9 and 5 × 10 -8 mol/m 2 /s.Solid black line give the median solution (D e = 5 × 10 -12 m 2 , DC 0 = 3.84 mol/m 3 , while dashed lines show extreme solutions.Minimum estimates of D e and DC 0 give maximum calculated reactions times (D e = 1 × 10 -12 m 2 , DC 0 = 3.34 mol/m 3 ), and maximum estimates of D e and DC 0 give minimum calculated reactions times (D e = 1 × 10 -11 m 2 , DC 0 = 4.34 mol/m 3 ).
FIgure 7. Plot showing the solution to the analytical one-component reactive transport equations.The duration it takes for the hematite reaction front to migrate 3.1 mm calculated from Equation5given the time for the reaction front to develop (t 0 ; Eq. 4).The shaded area represents the solutions for 427 < q < 896 m -1 , implying timescales of reaction between 4 and 20 yr and hematite dissolution rates between 1 × 10 -9 and 5 × 10 -8 mol/m 2 /s.Solid black line give the median solution (D e = 5 × 10 -12 m 2 , DC 0 = 3.84 mol/m 3 , while dashed lines show extreme solutions.Minimum estimates of D e and DC 0 give maximum calculated reactions times (D e = 1 × 10 -12 m 2 , DC 0 = 3.34 mol/m 3 ), and maximum estimates of D e and DC 0 give minimum calculated reactions times (D e = 1 × 10 -11 m 2 , DC 0 = 4.34 mol/m 3 ).

FIgure 8 .
FIgure 8. Diffusion-reaction modeling of gypsum-siltstone alteration zone in PHREEQC at 6 yr since start of injection; with (a) CO 2 -SO 4 -saturated infiltrating fluid (Model 1) and (b) CO 2 -poor SO 4 -saturated infiltrating fluid (Model 2).One thousand iterations of each model were run (gray lines) with the infiltrating fluid compositions randomly varied ±10% from the compositions presented in Table1.The black lines are the mean of all these models.The precipitation of illite in the CO 2 -saturated model matches the observed mineralogical profiles the best, suggesting the fluids responsible for the alteration were CO 2 -saturated.

FIgure 10 .
FIgure 10.Comparison between CO 2 -SO 4 -saturated PHREEQC model at 6 yr with an effective diffusivity (D e ) of 5 × 10 -12 m 2 and observed mineralogical profiles: (a) hematite profile; (b) dolomite profile; and (c) illite profile.Squares represent the average concentration over a 0.31 mm block across the alteration zone, these 26 points were determined from the QEMSCAN mineralogical maps.Model shown was the best fit to the observations of all simulation based on least squares.The modeling predicts the trends of the dissolution of hematite and dolomite upstream of the hematite reaction front, and growth of secondary illite adjacent to the gypsum bed.
measured the hydration waters of this gypsum bed at 158.48 m, the sample comprised 19.7% H 2 O. Figures3 and 4show 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: