Appropriate model design is critical for equilibrium calculations that predetermines the stability of mineral phases (Hassannayebi et al., 2019; Bo et al., 2021). To evaluate Fe-H2 interaction under such equilibrium conditions, we have selected water compositional data (Table 1) from a shallow brine system in the Murray Basin of southeastern Australia (Harrington et al., 2008). This was chosen to test different Fe speciation in a mildly acidic to circumneutral brine (Table 1) under near ambient conditions that mimic most of the natural brines of the world. As the pH range of such natural brines varies between ~ 4 to ~ 8, it is also necessary to constrain the Fe activity for any site before hydrogen injection, as described in the following sections. Temperature is set at a lower estimate of 50 \(℃\) to mimic physical conditions for most natural brines to approach the activation energy barriers for abiotic reactions. Although, as observed from the simulations, the system is little affected by a change in pressure; 1 MPa is set as the maximum estimate for our calculations (i.e. shallow system). Simulations were carried out using reductive dissolution for simple minerals like hematite, and magnetite under different initial Fe activities.
Activities of all the aqueous species were calculated using the Debye-Hückel formulation (B-dot) activity model within SpecE8 of Geochemist’s Workbench (GWB) v.17 (Bethke, 2022). Calculated model ionic strengths are approximately 1 M, the upper limit of the B-dot activity model (Xiong, 2006). The B-dot equation (Helgeson, 1969) has been developed to span a wide temperature range (up to 300°C). Thermodynamic calculations within GWB utilized the standard thermo database from Lawrence Livermore National Laboratory (LLNL) repository. Some major ion molalities are provided in Table 1 at the end of this section. In this study, we aim to delineate the general behavioral pattern of dissolved H2 in the presence of different concentrations and redox states of aqueous Fe; therefore, only a few representatives from major dissolved ion data seem sufficient. For reference, readers can access the entire database of the brine from the mentioned study (Harrington et al., 2008). Although present in trace quantities, Fe ions can show a relatively large variation in their concentrations in natural brines, ranging from micromolar to nanomolar scales (Trapp and Millero, 2007). Therefore, this is not considered a major cation species for the initial calculations. The effect of variation in Fe concentration also does not significantly alter the entire equilibrium state of the system, except for some typical Fe-bearing phases (e.g. hematite, magnetite, goethite, pyrite, and siderite). Dissolved Fe2+ concentrations are chosen to range within 1 milimolal (high endmember) to 1 micromolal (lower endmember) for quantification purposes (Supplementary Table 1), and Fe3+ is assumed to be solely governed by the dissolved hydrogen concentration. The oxygen fugacity (fO2) values were not readily available for the studied brine, although a wide range (-15 to -70) of values are reported from different regions worldwide. Such estimation is also supported by the fact that previous studies have reported Fe redox systematics can still be quite active in brine environments with very low dissolved oxygen concentrations, much less than 1 micromolal (Grenthe et al., 1992). In the cold, natural, and near-neutral/alkaline geothermal waters of Iceland, previous investigations have suggested that partial equilibrium between the Fe2+/Fe3+ were approached (Stefánsson et al., 2005). Due to the formation of several partially oxidized metastable Fe3+-Fe2+ aqueous complexes, fluctuations in Eh and pH occur during the hydrolysis of ferrous (Fe2+) and ferric (Fe3+) ions in aqueous solutions that governs the affinity for major thermodynamic changes in brines (Córdoba et al., 2016). Abiotic iron redox pairs can equilibrate rapidly within days in circumneutral pH conditions (St Clair et al., 2019), and for waters with pH\(\ge\)6, Fe approaches complete redox equilibrium (Kaasalainen et al., 2017). Formation of Fe(OH)2 and Fe(OH)+ species also controls the overall kinetics of Fe2+ to Fe3+ oxidation in near-neutral to mildly acidic pH ranges (~ 5–6). Similarly, soluble FeCl+ complexes also determine the redox kinetics and Fe2+ activities in natural brines (Millero, 1985; Trapp and Millero, 2007), where the latter is altered by anionic complexations. Aqueous silica concentration is set at 0.1 milimolal, a reasonable upper estimate for brines without disrupting the overall equilibria. The effect of silica on Fe-H2 systematics, with regard to the alteration of primary mineral phases, is studied elsewhere (Frost and Beard, 2007; Syverson et al., 2017) and is not the primary goal of this study. However, silica concentrations as high as 500 ppm is reported from many geothermal brines at surficial conditions (Gallup et al., 2003); thus, a maximum estimate of ~ 100 ppm is considered for understanding the behavior of the concerned species in the equilibrium studies. Activity diagrams are generated through thermodynamic phase relationships between different species and mineral stabilities are considered as forward precipitation reactions where the precipitating mineral is thermodynamically favored. In other words, the saturation index (SI) must be greater than 0 for precipitation to get initiated, which is calculated using-
SI =\(log \left(\frac{Q}{K}\right) \to \{<0 undersaturated\)
Q is the ion activity product, and K is the equilibrium constant for the relevant equilibrium reactions. This study focuses on pure mineral phases with the assumption of instantaneous precipitation/dissolution if the equilibrium is altered.
The model considers no microbial influence to reduce the number of variables and uncertainty. In one set of kinetic simulations, except for Fe, all other redox couples (bicarbonate-methane, nitrate-nitrogen) were disabled to test the Fe systematics in the absence of biology and at higher salinities for the aqueous media (Fanning, 2000; Truche et al., 2010). Fe redox and even Thermochemical Sulfate Reduction (TSR) are effective at low to near neutral pH (Kiyosu, 1980; Santana-Casiano et al., 2005; Zhang et al., 2008, 2012), and their Eh values often overlap. Because the onset and rate of TSR are influenced by several variables that differ from location to location, such as the composition of the available organic reactants, kinetic inhibitors and/or catalysts, anhydrite dissolution rates, wettability, as well as migration and diffusion rates of the major reactants toward one another, TSR does not have a clearly defined, minimum temperature (Goldstein and Aizenshtat, 1994; Machel, 2001), although it can occur at temperatures as low as ~ 25 oC (Mougin et al., 2007). Therefore, sulfate-bisulfide redox couple was decopuled in this study, representing the maximum estimate of such potential redox effects without microbial actions. This is also justified on the assumption that at near circumneutral pH, Fe3+-Fe2+ redox can get activated without TSR (Esnault et al., 2010; Kaasalainen et al., 2017). In static simulations, the extent of maximum reduction of dissolved sulfate can be seen from Fig. 2.
In another set of simulations where all the redox couples are in disequilibrium, relevant kinetic rate laws are applied to the Fe3+-Fe2+ redox pair where they are related through the redox reaction of EQ 3. At low temperatures, Fe3+ bearing species (hematite) was considered to dissolve reductively to release Fe2+ ions in the presence of gaseous H2. For this study, the abiotic reaction rate constant (at \({50 ℃, k}_{50}\)) at near acidic pH (~ 6) can be calculated through a simple Arrhenius-type relationship from the data obtained from Palandri and Kharaka (2004):
$$log\frac{{k}_{50}}{{k}_{25}}= -\frac{{E}_{a}}{R} \left(\frac{1}{348.15}- \frac{1}{298.15}\right)$$
\({E}_{a}=\) 66.2 KJ/mol, \({k}_{25}=\)4.07 \(\times\)10−10 mol/m2/s (for acid mechanisms)
In all cases, Fe oxides’ reductive dissolution depends upon partial H2(g) pressure to the first order (Fischer, 1987; Otake et al., 2010). This approach led to a kinetic rate constant ranging between 10− 9 to 10− 10 mol/m2/s, which is of the similar order of magnitude as reported in Torrent et al. (1987). Note that the results are not significantly altered when rate constants are varied within the mentioned range as the water-rock ratio is an important factor for the extent of such reactions.
Rate is calculated through React module of Geochemist’s Workbench (GWB) v.17 (Bethke, 2022) through:
$$r=A{k}_{50}\left(1-\frac{Q}{K}\right)$$
\(A=\) reactive surface area (taken as 20 cm2/g for hematite and magnetite)
\(Q\) = ion activity product and \(K=\)equilibrium constant
An initial dissolved CO2 was taken at 0.15 molal and considered temperature, pressure, and salinity (Yan et al., 2011). This choice is crucial as increasing CO2 activity would gradually stabilize carbonate phases like siderite. For the sake of simplicity in our kinetic simulations, effects of sulfate or bicarbonate ions are neglected.
Table 1
Basic chemistry model input parameters at T = 50\(℃\). Ion concentrations are taken from Harrington et al. (2008).
pH | Na+ (mol/kg) | Mg2+ (mol/kg) | K+ (mol/kg) | Ca2+ (mol/kg) | Cl−(mol/kg) | SO42−(mol/kg) | HCO3− (mmol/kg) |
5.9 | 0.542 | 0.058 | 0.006 | 0.008 | 0.652 | 0.039 | 0.88 |