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

Introduction

Storing anthropogenic carbon dioxide (CO₂) 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₂) 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₂ 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₂ stripped from the flue gas during the carbon capture process is likely to contain trace (<5%) co-containments, i.e., SO₂, H₂S, NOx, and O (Talman 2015). To purify the gas further is a costly process and hence, if CO₂ can be safely stored with these co-contaminants it will be economically beneficial. When SO₂ dissolves in groundwater it may form H₂SO₃, H₂SO₄, or H₂S depending upon redox conditions, and these compounds are stronger acids than dissolved carbonate species formed through CO₂ 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 CO₂ on cap rocks and reservoir rocks. These long exposure times may lead to porosity-permeability changes from mineral precipitation and dissolution. The study of natural analogs with S- and CO₂-bearing reservoir fluids are needed to constrain long-term predictions of reservoir and cap rock alteration during CO₂ storage.

Natural CO₂ and CO₂-SO₄-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₂. Scientific drilling adjacent to the Little Grand Wash fault in 2012 has allowed study of CO₂ transport mechanisms and long-term reactions between S-bearing, CO₂-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 CO₂- 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 CO₂-saturated at depth (~320 m) with minor S and CH₄ (Kampman et al. 2014).

Figure 1

Geological map of the Green River anticline showing locations of the Little Grand Wash and Salt Wash Graben normal fault systems, CO₂-springs and location of drill-hole CO2W55. Structure contours are for the top surface of the Navajo Sandstone, the main CO₂-bearing 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.

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₂-bearing fluids on clay-rich siltstones.

Regional geology

Near the town of Green River, natural CO₂-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 CO₂-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₂ 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 CO₂ leakage for at least 400 000 yr (Burnside et al. 2013). The intermittent travertine deposition is driven by the CO₂-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₂-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₂ 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., CO₂-S-bearing brines) must be more dense than the formation fluid. The modern fluids in the Entrada and Navajo Sandstone are CO₂- and SO₄-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 CHt (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 SO₄ contain ~0.5 mmol/L H₂S) coupled with the evidence for precipitation of pyrite in the basal portions of the Entrada and Carmel cap rocks. CO₂-bearing core experiments conducted by Purser et al. (2014) on red Entrada Sandstone showed H₂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³⁺ in the core sample to Fe²⁺ (mobilizing 15 mg/L of iron), thus bleaching the entire sample in the 6 month experiment. Whereas, using 5% CH₄ 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₂S rather than with CH₄ as the reducing agent. Given this evidence, H₂S is considered to be the primary reductant in the Green River CO₂ accumulation.

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

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) 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:

Table 1

Compositions (mol/L) of initial pore water (Turrero et al. 2006), CO₂-SO₄-saturated brine (measured sample DFS004 saturated with gypsum; Kampman et al. 2014), and CO₂-poor SO₄-saturated brine used in the PHREEQC reactive transport modeling

pH	Al3+	Ca2+	Mg2+	K+	Na+	SiO2	Cl-	SO42−	H2S	Alk	DIC
Initial pore water
7.52	0	13.6	1.63	0.43	23.9	0	30.4	23.9	0	0.76	0.8
CO2-SO4-saturated brine
5.1	0.01	24.1	10.1	7.2	112.5	0.1	84.9	20.7	0.5	64	884
CO2-poor SO4-saturated brine
5.1	0.01	24.1	10.1	7.2	112.5	0.1	84.9	20.7	0.5	2.5	2.5

1. Alk = alkalinity (mEq/L), DIC = dissolved inorganic carbon.

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

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.

Albite: NaAlSi₃O₈

Ankerite: CaFe_o.7Mg_0.3(CO₃)₂

Dolomite: CaMg(CO₃)₂

Hematite: Fe₂O₃

Illite: K_0.97(H₃O)_0.4Al_2.46Fe_0.3Mg_0.17Ti_0.03Si_3.07O_7.15(OH)₂·(H₂O)

Orthoclase/Microcline: K(AlSi₃O₈)

Quartz: SiO₂

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₂, SO₂, and H₂S that result in the removal of hematite, dolomite, orthoclase, and albite via stoichiometries such as:

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.

Hematite dissolution:

4Fe₂O₃ + H₂S + 14CO₂ + 6H₂O →8Fe²⁺ + SO42−CO3− Fe₂O₃ + 2HS^– + 4CO₂ + H₂O → Fe²⁺ + FeS₂ + 4HCO3−

Dolomite dissolution:

MgCa(CO₃)₂ + 2H⁺ → Mg²⁺ + Ca²⁺ + 2HCO3− MgCa(CO₃)₂ + 2H₂CO₃ → Mg²⁺ + Ca²⁺ + 4HCO3− MgCa(CO₃)₂ ↔ Mg²⁺ + Ca²⁺ + 2CO32−

Orthoclase/Microcline dissolution:

K(AlSi₃O₈) + 9H⁺ → K⁺ + Al³⁺ + 3H₃SiO₄

Albite dissolution:

NaAlSi₃O₈ + 9H⁺ → Na⁺ + Al³⁺ + 3H₃SiO₄.

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:

CO₂(g) + H₂S + albite + dolomite + orthoclase/microcline + hematite → illite + ankerite+Mg²⁺+Fe²⁺+Na⁺ + Ca²⁺ + HCO3− + SO42−

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², effective solute diffusivities, D_e, in the fluid phase on the order of 10^–11 to 10^–12 m² (Kampman et al. 2016), and a piezometric pressure gradient between 10⁴ and 10⁶ 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 (ω₀φ) of <6 × 10^–10 m/s for the siltstones, and Peclet numbers (Pe = ω₀φh/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³) with time t (s) and distance, x(m) can be written as:

ϕ∂C∂t=De∂2C∂x2+κfα(Ceq−C)(1)

where D_e is the effective diffusion coefficient, k_f is the mineral reaction rate (m/s), α is the mineral surface area (m²/m³), and C_eq is the solute concentration at equilibrium. The variation in the mineral molar concentration in the rock, φ_S(mol/m³) with time t, as controlled by surface dissolution, is modeled by a linear kinetics as:

∂φS∂t=−kfα(Ceq−C).(2)

Following Kampman et al. (2016), the following transformations to dimensionless variables are made:

x=lx′t=t′l2/DeC=Ceq−C′(Ceq−C0x=lx′(3)

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′∂t′=∂2C′∂x′2+NDC′(4)

∂φs′∂t′=−NDC′(5)

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:

ND=kfαl2De(6)

For a mildly reducing, S-bearing, CO₂-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²/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 ~10⁶ m²/m³. 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 (τ₀) it takes for hematite at the vein contact to become exhausted is described, assuming a constant hematite surface area, by:

τ0=φS∞VskfαΔC0(7)

where φS∞ is the dimensionless volume of hematite initially present. For times, t > τ₀, hematite is exhausted over distances x ≤ l.

Once this hematite dissolution front has formed (t > τ₀), the time it takes for the front to migrate an observed distance, l, is calculated by:

l+12ql2=1qτ0(t−τ0)(8)

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:

q=kfαDe.(9)

Note that the Damköhler number, N_D = (ql)².

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:

φs(x,l)=φS∞(1−e−q(x−1))(10)

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 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, D_e = 5 × 10^–12 m²/s, an estimated hematite surface area of 10⁶ m²/m³ 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²/s (gray shaded area; Fig. 7).

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 1σ St.dev. The initial hematite volume, φ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.

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 Equation 5 given the time for the reaction front to develop (τ₀; 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^–8mol/m²/s. Solid black line give the median solution (D_e = 5 × 10^–12 m², ΔC₀ = 3.84 mol/m³, while dashed lines show extreme solutions. Minimum estimates of D_e and ΔC₀ give maximum calculated reactions times (D_e = 1 × 10^–12 m², ΔC₀ = 3.34 mol/m³), and maximum estimates of D_e and ΔC₀ give minimum calculated reactions times (D_e = 1 × 10^–11 m², ΔC₀ = 4.34 mol/m³).

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.9 and the llnl.dat database (Parkhurst and Appelo 1999).

Key parameters and variables

The 3.1 mm wide gypsum-siltstone 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 (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.3Fe_0.7Ca(CO₃)₂] 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₂, pO₂, 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₂.

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₂-poor, SO₄-rich fluid likely relating the uplift of the Colorado Plateau after the Eocene, or (2) a CO₂-SO₄-saturated fluid relating to the recharge phase of the CO₂-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₂ content from saturated (Model 1) to CO₂-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₂-SO₄-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− H₂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 was varied and bounding fluid composition kept constant for the CO₂-saturated case only. Of the 1000 simulations in Models 1 and 2, 13 of the CO₂-saturated simulations and 3 of the 1000 CO₂-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₂S, and SO42− 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₂- SO42− -saturated infiltrating fluid. Critically the CO₂-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₂-SO42−-saturated model, suggesting the alteration is associated with the CO₂-accumulation. Additionally, the velocity of the hematite reaction front is directly related to the concentration of H₂S and the H₂S/SO42− ratio in the infiltrating fluid. Faster front velocities are associated with a greater H₂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 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₂ 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²; 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² 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 theinference of timescales.

Figure 8

Diffusion-reaction modeling of gypsum-siltstone alteration zone in PHREEQC at 6 yr since start of injection; with (a) CO₂-SO₄-saturated infiltrating fluid (Model 1) and (b) CO₂-poor SO₄-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 Table 1. The black lines are the mean of all these models. The precipitation of illite in the CO₂-saturated model matches the observed mineralogical profiles the best, suggesting the fluids responsible for the alteration were CO₂-saturated.

Figure 9

Diffusion-reaction modeling of gypsum-siltstone alteration zone in PHREEQC at with a CO₂-SO₄-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², 5 × 10^–12 m², and 1 × 10^–11 m², respectively. The plot shows the effect diffusivity has on the propagation speed of the reaction fronts.

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²/s. A comparison at a 6 yr time step between the CO₂-SO₄-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₂-poor fluids. Discrepancies between the observed mineral profiles and the CO₂-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₂-SO₄-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.

Figure 10

Comparison between CO₂-SO₄-saturated PHREEQC model at 6 yr with an effective diffusivity (D_e) of 5 × 10^–12 m² 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.

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, 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₂-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^–12m^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, 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₂-S fluid-rock reactions will be.

Other studies of CO₂ 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₂ alteration resulted in the formation of Fe²⁺- and K⁺-enriched illites. Overall this study agrees with those of Busch et al. (2014) and Kampman et al. (2016) who found that CO₂-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₂ 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₂-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.