Geofluids
Volume 2018 (2018), Article ID 9475741, 20 pages
https://doi.org/10.1155/2018/9475741

Research Article

Evaluation of the Potential for Dissolved Oxygen Ingress into Deep Sedimentary Basins during a Glaciation Event

Sergio A. Bea,¹ Danyang Su,² K. Ulrich Mayer,² and Kerry T. B. MacQuarrie³

¹CONICET-IHLLA, República de Italia 780, C.C. 47, Azul, Buenos Aires B7300, Argentina
²Department of Earth, Ocean and Atmospheric Sciences, University of British Columbia, 2207 Main Mall, Vancouver, BC, Canada V5T 1Z4
³Department of Civil Engineering, University of New Brunswick, P.O. Box 4400, Fredericton, NB, Canada E3B 5A3

Correspondence should be addressed to Sergio A. Bea; sabea@faa.unicen.edu.ar

Received 24 August 2017; Revised 9 February 2018; Accepted 18 February 2018; Published 24 April 2018

Academic Editor: Tianchyi Yeh

Copyright © 2018 Sergio A. Bea et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Geochemical conditions in intracratonic sedimentary basins are currently reducing, even at relatively shallow depths. However, during glaciation-deglaciation events, glacial meltwater production may result in enhanced recharge (Bea et al., 2011; and Bea et al., 2016) potentially having high concentrations of dissolved oxygen (O₂). In this study, the reactive transport code Par-MIN3P-THCm was used to perform an informed, illustrative set of simulations assessing the depth of penetration of low salinity, O₂-rich, subglacial recharge. Simulation results indicate that the large-scale basin hydrostratigraphy, in combination with the presence of dense brines at depth, results in low groundwater velocities during glacial meltwater infiltration, restricting the vertical ingress of dilute recharge waters. Furthermore, several geochemical attenuation mechanisms exist for O₂, which is consumed by reactions with reduced mineral phases and solid organic matter (SOM). The modeling showed that effective oxidative mineral dissolution rates and SOM oxidation rates between 5 × 10⁻¹⁵ and 6 × 10⁻¹³ mol dm⁻³ bulk s⁻¹ were sufficient to restrict the depth of O₂ ingress to less than 200 m. These effective rates are low and thus conservative, in comparison to rates reported in the literature. Additional simulations with more realistic, yet still conservative, parameters reaffirm the limited ability for O₂ to penetrate into sedimentary basin rocks during a glaciation-deglaciation event.

1. Introduction

In intracratonic sedimentary basins, reducing geochemical conditions are widespread and commonly dominate even at shallow depths close to the ground surface. The maintenance of reducing conditions is often cited as one of the desirable features for deep geological repositories (e.g., [1, 2]). In general, oxygen (O₂) ingress to repository depth is considered an exclusion criterion for the construction of deep geologic repositories in view of the potential for oxidative dissolution reactions and the corrosion of copper canisters containing the spent fuel (e.g., [3]). One of the explanations for the rapid development of reducing conditions with increasing depth is that O₂ in the recharge water is consumed due to interactions with redox-buffering minerals (e.g., pyrite and chlorite) present as coatings on the aquifer matrix and solid organic matter (SOM) [2, 3]. However, over long time periods, glaciation and deglaciation events may affect groundwater flow patterns and the distribution and rate of recharge. For instance, during deglaciation events, melting of ice sheets may result in enhanced recharge of glacial meltwater (e.g., [4–9]). McIntosh et al. [7] summarized geochemical and isotopic evidence for freshwater ingress in North American Paleozoic basins as a consequence of glaciation cycles. For example, McIntosh et al. [7] and references therein showed that freshwater significantly diluted saline brines in the Upper Devonian Antrim Shale across the northern margin of the Michigan Basin as evidenced by a decrease (up to 80%) in bromide concentrations to a depth of 300 m [6]. The source of the freshwater recharge was inferred to be subglacial meltwater, because the ¹⁴C ages correspond to periods of ice cover, and the low δ¹⁸O values provide evidence for mixing of ice sheet derived waters with saline brines and modern precipitation [10, 11]. Several other associated studies reported evidence for recharge and storage of glacial meltwaters in the Michigan Basin, documented by low salinity, ¹⁴C ages ranging from modern to >50 ka, and oxygen isotope values ranging from −18 to −13 [4, 12–14]. In addition, noble gas signatures have also been used to infer the ingress of cold glacial meltwaters [15, 16].

These recharge waters may also be enriched with dissolved O₂ relative to present-day interglacial conditions by a factor of 3–5 (e.g., [2, 20]), providing an enhanced oxidation capacity. Oxygen enrichment can be attributed to air entrapment in ice during periods of the glacial advance. During melting of an ice sheet entrapped air may be released into meltwater under elevated hydraulic head conditions, explaining the enriched dissolved O₂-concentrations [2, 20]. It is therefore of interest to evaluate the probability for enhanced dissolved O₂ ingress and migration under these changed recharge conditions. This information can be used to improve the formulation of conceptual and quantitative models for geological repository safety assessment.

The key parameters controlling O₂ ingress and attenuation in crystalline rocks have been previously identified and evaluated using reactive transport modeling (e.g., [2, 3, 20]). In crystalline rocks, the penetration depth of O₂ depends largely on the spatial distribution and interconnectivity of subsurface pathways within the fractured porous medium, controlling advective flow. Sidborn and Neretnieks [3] determined that O₂ may travel long distances in the absence of readily degradable organic matter but also acknowledged that large uncertainties regarding key model parameters remain, leading to the conclusion that the results must be treated with caution pending more accurate and validated data. Although these studies provide information on the redox stability of crystalline environments, we are unaware of observational or modeling studies focused on O₂ ingress, migration, and attenuation in deep sedimentary basins, in particular over long time frames involving glaciation-deglaciation cycles. Deep sedimentary basins are complex systems requiring the use of a thermo-hydro-chemical and mechanical approach to simulate regional groundwater flow and transport processes in a more realistic way (e.g., see [8, 9]).

The objective of this paper is to evaluate dissolved O₂ ingress and attenuation in a hypothetical sedimentary basin during a glaciation/deglaciation event by exploring key physical and geochemical parameters that affect the persistence of O₂, thereby contributing to an improved understanding of the phenomena governing redox stability of deep geological repositories in sedimentary environments. To achieve this objective, Par-MIN3P-THCm [8, 9, 21–23] is used to perform reactive transport modeling on the basin scale combined with a parameter sensitivity analysis. Transient surface boundary conditions are imposed on the upper part of the model domain to mimic ice sheet advance and retreat (e.g., [9, 24]). Processes included in the current model are nonisothermal density-driven flow and transport, a simplified model accounting for vertical mechanical deformation, as well as geochemical reactions (e.g., aqueous complexation and mineral dissolution and precipitation including the interaction of dissolved O₂ with redox-buffering solid phases such as pyrite, chlorite, biotite, and SOM). We acknowledge that the selected mineral phases and SOM represent only a small subset of reduced phases present in sedimentary rocks. However, they represent some of the most abundant and redox active phases, providing an adequate framework to assess O₂ ingress and attenuation. Regrettably, it is not possible to directly determine oxidative dissolution and SOM oxidation reaction rates to allow for deterministic simulations. Rates are affected by mineral abundance, surface normalized rate constants, effective reactive surface areas, are often microbially mediated, and differ as a function of lithology. In the absence of quantitative knowledge of O₂ consumption rates in sedimentary basins, the model is used to identify the effective reaction rates required to restrict O₂ ingress to depths of less than 200 m over the duration of an entire glaciation-deglaciation cycle and to compare these rates to published data on oxidative mineral dissolution rates and oxidation of SOM. This approach allows to evaluate the effectiveness of the redox-buffering capacity of sedimentary basins during periods of glacial meltwater ingress. A relatively shallow depth of 200 m was chosen as tolerable for O₂ ingress because deep geologic repositories in sedimentary basins are being planned at significantly greater depths (e.g., 500 to 700 m), and 200 m thus provides a conservative redox buffer depth relevant to such repositories. In addition, the sensitivity of O₂ ingress towards effective rates and O₂ content in recharge water is investigated.

2. General Features of Deep Sedimentary Basins and Impact of Glaciation Cycles on the Groundwater System

Intracratonic sedimentary basins in North America and Eurasia are commonly characterized by Paleozoic sequences of limestones/dolostones interbedded with shales and sandstones and in some cases evaporites (e.g., [1, 4–7]). In the Illinois, Michigan, and Appalachian basins in North America, for instance, the basal Cambrian-Ordovician aquifers are principally comprised of sandstones and carbonates (dolomite and limestone) and confined by shales (e.g., [4]), whereas the Silurian-Devonian aquifers are primarily composed of dolostones and limestones confined by the overlying Upper Devonian black shales, Mississippian grey shales, and siltstones. Furthermore, Silurian-Devonian aquifers in the Michigan and Appalachian basins contain thick sequences of interbedded halite [NaCl] and anhydrite [CaSO₄], which are nearly impermeable (e.g., [4]). Formational waters in these deep Paleozoic sedimentary basins are brines (e.g., [4, 6, 7]) representing in some cases (e.g., in the Michigan Basin) the most saline fluids in the Earth’s crust [5, 6]. Notably, brines within these basins are of different types (CaCl₂ and NaCl-type), showing a spatially variable density distribution [6, 8, 9].

The formation and melting of large ice sheets have been proposed as a mechanism for influencing regional-scale groundwater flow systems during the Late Pleistocene in sedimentary basins located in the interior of the North American and Eurasian cratons. Hydrochemical and isotopic evidence suggests that glacial meltwaters may have displaced high salinity fluids and thus stored freshwater resources on the continents [4, 6, 7, 24–26]. During ice sheet loading, small changes in pore volume as the rock expands or compresses may result in substantial localized changes in the pore fluid pressure [24]. Glacial loading and unloading are considerably faster processes than sedimentation and erosion [24], while these changes in fluid pressure quickly dissipate in highly permeable units, the effects of loading and unloading cycles could induce anomalous fluid pressures that require thousands of years to return to equilibrium conditions in low-permeability units like shales (e.g., [8, 24, 27]). These effects on the flow regime may also affect the ingress and fate of dissolved O₂ into the sedimentary units.

Under present-day conditions, advective solute migration in shallow permeable aquifer units, including transport of dissolved O₂ contained in the surficial recharge water, is limited to shallow depths [8, 9]. In addition, dissolved O₂ could potentially be consumed by reactions with reduced mineral phases (e.g., pyrite, chlorite, or biotite) and SOM, further contributing to O₂-depleted conditions near the ground surface. The depth to which O₂ can migrate will depend on the interplay between altered fluid circulation patterns, the O₂ concentration in recharge water, and the abundance and reactivity of SOM, as well as sulfide and Fe(II)-bearing mineral phases in the sedimentary units.

3. Simulation Approach and Material Property Assignment

Reactive transport simulations were carried out with the code Par-MIN3P-THCm [8, 9, 22, 23]. This code includes the modeling capabilities of previous versions (i.e., MIN3P [21]; MIN3P-D, [28]) but has been enhanced to account for the coupling of thermo-hydro-chemical and mechanical processes [24]. Specific enhancements include calculation of ion activity corrections in high ionic strength solutions using the Harvie-Möller-Weare (HMW) model [29] based on the Pitzer [30] equations; fluid density calculations based on fluid volumetric predictions derived from the Pitzer equations according to Monnin [31]; one-dimensional hydromechanical coupling [27] due to ice sheet advance and retreat; and a numerical formulation to couple energy (heat) transport (following Voss and Provost; [32]) with fluid and solute transport (i.e., density-driven flow controlled by thermohaline convection). In addition, the code has been parallelized to increase its computational efficiency and reduce execution times [23].

3.1. Conceptual Model and Physical Parameters

The conceptual model and parameterization for the hypothetical sedimentary basin considered here include the key features of several of the main Paleozoic sedimentary basins in North America. Physical and geochemical model parameters were mainly derived from previous studies in the Michigan and Illinois basins as described in detail in Bea et al. [8, 9] (i.e., porosity, horizontal and vertical hydraulic conductivity, elastic modulus, Poisson’s ratio, volumetric heat capacity, and bulk density).

The simulated basin domain (Figure 1) extends over a horizontal distance of 400 km, has a maximum depth of 4 km, and is composed of a sequence of carbonates (dolostones and limestones, Dol1, Dol2, Dol3, and Lim1), interbedded by sandstones (Sand1, Sand2, Sand3, and Sand4, constituting the main aquifers), and shales (Sh1, Sh2, and Sh3, constituting the main confining units). These units overlay a crystalline basement (Cr). A weathered zone in the crystalline basement (Crw) is assumed to be in direct contact with the sedimentary units. Interbedded evaporites (Ev) are also included within the dolostone units. Note that thin layers (e.g., <50 m) of unconsolidated glacial drift are not included in the current conceptual model. These units are removed and reform during periods of glacial advance and retreat and are not relevant for considering meltwater and solute ingress into the deeper sedimentary rock formations.

Figure 1: Domain, geology, and stratigraphy used to simulate the impact of continental glaciation on O₂ ingress and attenuation in a hypothetical sedimentary basin.

The main physical parameters (e.g., porosity and hydraulic conductivity) are considered to decrease with depth following an exponential relationship based on the approach of Bahr et al. [33]. The depth dependence of the hydraulic conductivities was captured by applying the Carman-Kozeny relationship using the depth-dependent porosities and a hydraulic conductivity determined at a reference depth. The depth dependence of other parameters, such as the specific storage coefficient and mechanical loading parameters, was determined in a similar fashion. This approach captures the general trends in parameter variation as a function of depth and is well suited for conceptual simulations (for more details, see [9]).

3.2. Geochemical Parameters and Initial Conditions

Formation waters present in Paleozoic sedimentary basins in North America possess chemical compositions that vary between meteoric water and evaporated seawater [4, 6, 19, 34, 35]. Chemical components included in the present reactive transport simulations are H, Ca, Mg, Cl, SO₄, K, Total Inorganic Carbon (TIC), Na, Br, O₂, and H₂O. For the sake of simplicity, the components Al, Fe(II), Fe(III), and H₄SiO₄ are only considered as part of mineral phases but are neglected for solute transport, since they do not significantly affect solution density and aqueous speciation. The complete set of homogeneous and heterogeneous geochemical reactions included in the simulations is shown in Table 1. Thermodynamic constants for the reactions were taken from Wolery and Daveler [17], whereas the selectivity coefficients for cation exchange reactions were obtained from Appelo and Postma [18]. The virial coefficient database used for the HMW model [29] (at 25°C and 1 atm) was originally developed as a part of the Yucca Mountain project [36]. The fluid density calculations were also based on Pitzer’s equations [31, 37] for consistency with the thermodynamic model.

Table 1: Geochemical reactions considered in the illustrative sedimentary basin example.

Ref.

Homogeneous reactions

Cation-exchange reactions

Dissolution/precipitation reactions
Enthalpy
kcal mol⁻¹

Equilibrium constant taken from Wolery and Daveler [17]. Selectivity coefficients taken Appelo and Postma [18]. Not controlled by equilibrium [2].

Mineral phases considered include halite [NaCl], anhydrite [CaSO₄], calcite [CaCO₃], and dolomite [CaMg(CO₃)₂], together with the redox-buffering solids pyrite [FeS₂], biotite [K(Mg₂Fe)(Si₃Al)O₁₀(OH)₂], chlorite [Mg₂Fe₃Al₂Si₃O₁₀(OH)₈], and SOM (represented stoichiometrically as CH₂O). These redox-buffering solids were chosen based on their typical abundance and level of reactivity. Pyrite reacts readily with O₂ (e.g., [38]), while chlorite and biotite also belong to the group of the most reactive Fe-bearing silicate minerals [2]. Organic matter is commonly present in sedimentary rocks, considering the depositional environment. Although this approach provides a highly simplified representation of rock composition in sedimentary basins, it allows making conservative assumptions regarding the abundance of reduced phases in these rocks, which is adequate for this study. Initial mineral contents for the different units are tabulated in Table 2. The initial spatial distribution of volumetric fractions for calcite, dolomite, halite, and anhydrite is identical to those shown in Bea et al. [9], whereas the distribution of the redox-buffering solid phase is shown in Figure 2. Note that calcite, dolomite, and SOM are considered as primary phases in the entire domain, except for the evaporite units (Ev, Figure 1). In the evaporite units, only halite and anhydrite are considered as primary phases. In addition, pyrite is considered to be initially present in shale and sandstone units (Figure 2(a)), while biotite and chlorite are present in the crystalline basement and in the weathered layer (Figure 2(c)). Chlorite is also assumed to be present in the shale units (Figure 2(d)).

Table 2: Initial mineral contents (%) considered in the reactive transport simulations (for reference see Bea et al., [9]).

Mineral
Cr
Crw
Sand1
Lim1
Sand2
Sh1
Dol1
Ev
Sh2
Sand3
Dol2
Sand4
Sh3
Dol3

Calcite
0
0
9
89
9
9
10
4
9
9
10
9
9
10

Dolomite
0
0
1
10
1
1
89.99
36
1
1
89.99
1
1
89.99

Halite
0
0
0
0
0
0
0
50
0
0
0
0
0
0

Anhydrite
0
0
0
0
0
0
0
10
0
0
0
0
0
0

Pyrite
0
0
0.1
0
0.1
2
0
0
2
0.1
0
0.1
0
0

Biotite
4
4
0
0
0
0
0
0
0
0
0
0
0
0

Chlorite
0.7
0.7
0
0
0
10
0
0
10
0
0
0
0
0

Organic matter
0
0
0.1
1
0.1
0.4
0.1
0
0.4
0.1
0.1
0.1
10
0.1

Nonreactive minerals
95.3
90
90.0
0
90
90
0
0
90
90
0
90
90
0

Based on rock composition for Queenston, Georgian Bay, and Blue Mountain Formations reported by Mazurek [1].

(a) Pyrite

(b) Solid Organic Matter (SOM)

(c) Biotite

(d) Chlorite

Figure 2: Initial distribution of the volumetric contents of redox-buffering minerals. (a), (b), (c), and (d) are pyrite, organic matter, biotite, and chlorite contents, respectively. Vertical exaggeration 100 : 1.

The aforementioned redox-buffering mineral phases (i.e., pyrite, biotite, and chlorite) and SOM were modeled as irreversible kinetic reactions in the presence of O₂. The rate expression for the oxidative dissolution of pyrite is derived from the work of Weaver et al. [12], assuming that only the O₂ pathway is active, with a square root dependence on O₂. Considering the limited dependence on pH and relatively stable pH conditions in sedimentary basins, the rate expression was simplified by evaluating the rate constant for a pH = 7 to yield

where [mol dm⁻³ bulk s⁻¹] is the oxidative dissolution rate of pyrite, is the dissolved oxygen concentration [mol l⁻¹], is the kinetic rate constant evaluated at pH = 7 [mol0.5 l0.5 m⁻² mineral s⁻¹], is the mineral volumetric fraction of pyrite [dm³ mineral dm⁻³ bulk], is the density of pyrite [g dm⁻³ mineral], and is reactive surface area of pyrite [m² g⁻¹].

The kinetic rate expressions for chlorite and biotite have been adopted from reactive transport simulations in crystalline rock environments ([2] and references therein) and are described by the rate expression:

where [mol dm⁻³ bulk s⁻¹] is the oxidative dissolution rate of chlorite or biotite, is the corresponding kinetic rate constant [mol m⁻² mineral s⁻¹], is the mineral volumetric fraction of the phase considered [dm³ mineral dm⁻³ bulk], is the relevant mineral density [g dm⁻³ mineral], and is the reactive mineral surface area [m² g⁻¹].

The oxidation of SOM is represented by an expression similar to (2), using an effective rate coefficient, but also depending on the abundance of SOM to account for the distribution of organic matter in the different sedimentary units:

where [mol dm⁻³ bulk s⁻¹] is the oxidation rate of SOM, is the effective rate coefficient [mol dm⁻³ SOM s⁻¹], and is the volumetric fraction of organic matter present in the solid phase [dm³ SOM dm⁻³ bulk].

The kinetic parameters used in the simulations are summarized in Table 3. Using the volume fractions listed in Table 2, this approach accounts for the presence and abundance of reduced mineral phases and SOM in the various sedimentary units, when quantifying O₂ consumption rates. For the base case simulation, the reactive surface areas for pyrite, chlorite, and biotite and the effective rate coefficient for SOM were set to achieve oxidative dissolution rates that limit dissolved O₂ ingress with concentrations above 1 × 10⁻⁹ mol l⁻¹ (corresponding to 3.2 × 10⁻⁵ mg l⁻¹) to a depth of 200 m. As indicated above, this depth was selected to conservatively define a tolerable depth of O₂ ingress with regard to deep geologic repositories.

Table 3: Kinetic parameters and solid properties for the redox-buffering minerals used in base case simulation.

Redox-buffering mineral

mol l⁻¹
Density
g dm³ mineral
Reactive surface area
m² g⁻¹

Pyrite (diss.)
10⁻⁹
-
5,020
⁽⁶⁾5 × 10⁻⁵

Biotite (diss.)
⁽⁴⁾10−12.25
10⁻⁷
2,930
⁽⁶⁾5 × 10⁻⁵

Chlorite (diss.)
⁽⁴⁾10−12.25
10⁻⁷
2,600
⁽⁶⁾5 × 10⁻⁵

Organic matter (diss.)
(5,6)5 × 10⁻¹²
5 × 10⁻⁵
-
-

derived from Weaver et al. [12]; following Spiessl et al. [2] and references therein; in [mol0.5 l0.5 m⁻² minerals s⁻¹]; in [mol m⁻² mineral s⁻¹]; in [mol dm⁻³ SOM s⁻¹]; set for base case.

The model was run for 30,000 years using present-day boundary conditions to obtain a quasi-steady-state representation of conditions during an interglacial period. Pressure head, Darcy velocity, and fluid density distributions obtained during this simulation were used as the initial conditions for the subsequent glaciation/deglaciation simulation (following the approach of Bea et al. [9]). The initial temperature distribution at the top of the basin linearly varies from 0 to 10°C from the right to left side of the basin. This temperature distribution is representative of the north-south geographic temperature variation during the interglacial period, if the basin is located in northern latitudes. Because Par-MIN3P-THCm uses the Van’t Hoff and Arrhenius equations to account for the thermal dependence of the equilibrium constants and kinetic rates [8, 9, 21], respectively, the equilibrium constants were locally adjusted in each cell of the mesh to reflect the initial temperature distribution. Initial chemical water compositions were computed in agreement with these conditions assuming thermodynamic equilibrium with the primary mineral phases calcite and dolomite, in addition to gypsum and halite in the units with evaporites.

3.3. Discretization and Boundary Conditions for the Glaciation Scenario

The computational domain was discretized using 45,000 cells; 450 cells evenly distributed in the (horizontal) direction, and 100 cells distributed in the (vertical) direction. The upper part of the mesh was refined in order to better resolve the processes occurring in the shallow part of the aquifers.

The glaciation scenario boundary conditions for flow, energy, and solute transport are the same as used by Bea et al. [9], and they are shown in Figure 3. Note that, for flow, impervious boundaries were assumed for the two sides and the bottom of the domain. In order to illustrate the potential hydrogeological and geochemical perturbations generated by ice sheet movement over the sedimentary basin described above, a simplified glaciation scenario was applied along the top boundary (Figure 3(a)). The glaciation scenario consists of four stages that generically represent the last glacial maximum episode of ice advance and retreat in North America [39], following an approach similar to that used for previous illustrative basin scale simulations (e.g., [24]): Stage I corresponds to a phase of ice sheet accumulation (i.e., loading, for a period of 12,500 years), Stage II represents a constant ice sheet thickness for a period of 5,000 years, Stage III corresponds to the phase of ice sheet melting and retreat (i.e., unloading, for a period of 5,000 years), and Stage IV is representative of conditions with no ice present (i.e., return to present-day surface boundary conditions). Note that the loading/unloading cycle shown in Figure 3(a) is asymmetric. It is based on paleoclimatic reconstructions using isotopic data as well as ice sheet models that usually suggest a more rapid ice sheet retreat (driven by rising temperatures) than its previous accumulation (mostly driven by precipitation rates) (e.g., [24, 39]). Taking this asymmetric evolution into account, anomalous fluid pressures in recently glaciated areas should be expected, including regions of overpressure and underpressure (e.g., see [24]). In the current glaciation scenario, the maximum ice sheet thickness and extent were 2,000 m and 440 km, respectively, which implies complete ice sheet coverage of the basin during the Stage II glacial maximum.

(a) Evolution of maximum ice sheet thickness and extension

(b) Boundary conditions, Stages I and II

(c) Boundary conditions, Stage III

(d) Boundary conditions, Stage IV

Figure 3: Reproduced from [9] with permission. (a) Glaciation scenario used to specify boundary conditions along the top boundary of the model domain. (b), (c), and (d) Model boundary conditions for flow, solute, and energy transport equations. h: hydraulic head [m]; : temperature [°C]; : concentration [mol L⁻³]; : volumetric recharge [L⁻³ T⁻¹]; : solute mass flux [mol T⁻¹]; and : heat flux per unit of surface [E T⁻¹ L⁻²].

At the beginning of the simulation (and also for Stage IV, Figure 3(d)), a prescribed head representing topographically driven flow was imposed on the upper boundary. In this case, solute mass fluxes were based on meteoric water composition and heat fluxes were calculated based on a prescribed temperature distribution that was assumed for the upper boundary of the domain (Figure 3(d)). A constant heat flux of 0.05 W m⁻² was imposed at the bottom of the domain (e.g., [24]) (see in Figure 3) during all stages of the simulations.

During ice accumulation (i.e., Stage I), prescribed heads representing topographically driven flow were imposed on the periglacial area, the extension of which depended on the glaciation stage (Figure 3(a)), whereas a cold-based (nonmelting) condition, represented by a zero-flux boundary condition, was assumed under the ice sheet. The latter assumption was based on paleoclimate reconstructions in ice-covered regions during the last glaciation, which implied that groundwater recharge was prevented by overlying glaciers [24] (Figure 3(b)). A prescribed temperature of 0°C was imposed under the ice sheet. Note that, during this initial stage, the ice sheet progressively advanced over the basin from right to left, thus inducing a mechanical load on the underlying sedimentary units (Figure 3(a)); by 12,500 years the basin was assumed to be entirely glaciated.

During Stage III (Figure 3(a)), warm-based (i.e., melting) conditions were assumed beneath the retreating ice sheet. Under such conditions, a direct hydraulic connection between the glacial meltwater and the land surface is established (Figure 3(c)). A fixed hydraulic head corresponding to 95% of the ice sheet thickness was therefore specified [9, 24] (Figure 3(a)). Boundary solute fluxes beneath the ice sheet during Stage III were computed based on the water flux and an oxygen-enriched meltwater chemical composition (i.e., Solution 5 in Table 4).

Table 4: Chemical composition of brines and brackish waters considered in the reactive transport simulations (chemical compositions taken from Hobbs et al. [19]).

Formational and boundary water compositions

(4)
(5)

pH
6
6
6
6
6

⁽⁶⁾Ca
1.5
0.77
0.2
0.19
6.2 × 10⁻⁴

⁽⁶⁾Mg
6.5 × 10⁻³
4 × 10⁻²
8.6 × 10⁻³
10⁻²
8.2 × 10⁻⁴

⁽⁶⁾Cl
6
4.15
6.24
3.61
2.9 × 10⁻³

⁽⁶⁾SO₄
5.4 × 10⁻⁴
4 × 10⁻³
1.1 × 10⁻²
5.2 × 10⁻⁵
2.6 × 10⁻⁴

⁽⁶⁾K
3.4 × 10⁻²
5.7 × 10⁻²
0.07
0.12
6.4 × 10⁻⁵

(6,7)TIC
1.2 × 10⁻⁴
3.2 × 10⁻⁴
5.2 × 10⁻⁴
1.8 × 10⁻³
1.7 × 10⁻⁵

⁽⁶⁾Na
2.85
2.49
5.77
3
5.2 × 10⁻⁴

⁽⁶⁾Br

9.7 × 10⁻³
7.4 × 10⁻³
5.6 × 10⁻³
4 × 10⁻⁵

⁽¹⁰⁾O₂(aq)
-
-
-
-
1.0 (or 0.21)

⁽⁶⁾Fe(III)
-
-
-
-
-

⁽⁶⁾Al(III)
-
-
-
-
-

Chemical signature
Ca-Na-Cl
Na-Cl
Na-Cl
Na-Cl
-

⁽⁶⁾Ionic strength
7.6
4.97
6.48
3.82
5.2 × 10⁻³

⁽⁸⁾TDS
0.34
0.24
0.37
0.21
2 × 10⁻⁴

⁽⁹⁾Density
1217
1150
1204
1128
997

Units G, Gw, Sand1, and Sand2; units Lim1, Sh1, and Dol1; unit Ev; ⁽⁴⁾units Sh2, Sand3, Sand4, Sh3, Dol2, and Dol3; ⁽⁵⁾meteoric water; ⁽⁶⁾units in [mol l⁻¹]; ⁽⁷⁾total Inorganic Carbon; ⁽⁸⁾units in [kg l⁻¹]; ⁽⁹⁾units in [kg m⁻³]; ⁽¹⁰⁾partial pressure of O₂(g) in atm.

The conceptual model for glaciation-deglaciation implies that recharge of oxygenated water occurs: (a) during interglacial periods, (b) during the period of glacial advance in unglaciated regions in front of the advancing ice sheet, and (c) during glacial retreat in unglaciated regions and beneath the retreating ice sheet. During Stages I and II, corresponding to the period of glacial advance and glacial maximum, respectively, no O₂ ingress occurs beneath the ice sheet because cold-based conditions were assumed. During glacial retreat, elevated O₂ concentrations were assumed in the meltwater beneath the ice sheet.

4. Evaluation of O₂ Ingress and Attenuation

To investigate the impact of a glaciation cycle and uncertain geochemical parameters on O₂ ingress and consumption, a base case and three sensitivity scenarios were simulated. The main features of these scenarios are summarized in Table 5. The relevant rate parameters for the base case simulation were set such that O₂ ingress was limited to 200 m at all locations during the entire glaciation-deglaciation cycle for reasons discussed above. Sensitivity analyses were then used to evaluate O₂ ingress for slightly less conservative and more realistic assumptions. The results of the base case scenario are described in detail, while the sensitivity scenarios are used to provide additional information on the effect of O₂ consumption rates and O₂ concentrations in source glacial meltwater on depth of O₂ penetration.

Table 5: Description of the simulation scenarios included in the sensitivity analysis.

Scenario Name
Rate parameters
O₂ contents in meltwater

Base case
, from Table 3
1

Scenario 1
, increased 10-fold
1

Scenario 2
, from Table 3
0.21

Scenario 3
, increased 10-fold
0.21

Partial pressure of O₂(g) in atm corresponds to 60 mg l⁻¹ for 1 atm and 12.6 mg l⁻¹ for 0.21 atm.

5. Results and Discussion of the Base Case Simulation

For the base case scenario, the evolution of flow, thermal and hydrogeochemical conditions (excluding redox chemistry) are described in detail by Bea et al. [9]. However, to provide the necessary context, we summarize here the key aspects of these modeling results and simultaneously enhance the focus by including results related to the O₂ attenuation and redox stability of the basin.

5.1. Fluid Pressure, Flow, Solution Density, and Thermal Evolution

As a consequence of dense brines at depth and the near horizontal layering of the hydrogeological units that form the basin (note the significant vertical exaggeration, for instance, in Figure 4) the active flow regime is generally restricted to shallow regions (less than 300–400 mbgs) for the full duration of the simulation. In this context, the region affected by active flow is related to the depth of penetration of low salinity waters.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4: Simulation results for the illustrative sedimentary basin simulation. (a) and (b) Spatial distribution of point pressure heads at 20,000 years, and its temporal evolution at observation points located in a sandstone aquifer (Sand2, observation point located at 2,800 m depth) and in a shale aquitard (Sh1, observation point located at 2,300 m depth). (c) and (d) Spatial distribution of density at 20,000 years, and its temporal evolution at observation point located in sandstone aquifer (Sand4 observation point located at 300 m depth). (e) and (f) Spatial distribution of temperature at 20,000 years, and its temporal evolution at observation points located in sandstone aquifer (Sand1, depths indicated on figure).

Point pressure heads in the most deformable units increase due to ice sheet loading (e.g., Stage I, ice accumulation) throughout the basin. The magnitudes of the pressure head changes depend primarily on the hydraulic conductivity and the hydromechanical loading efficiency coefficient assigned to the various hydrogeological units (Figures 4(a) and 4(b)). Pressure heads increase beneath the retreating ice sheet in shallow permeable units (see results at 20,000 years, Figure 4(a)) as a consequence of the hydraulic connection between the now warm-based melting ice sheet and these sedimentary units (see Stage III in Figure 4(b)). Some of the highest Darcy velocities are generated in the shallow aquifers during the early stages of glacial retreat (e.g., on the order of 10⁻¹ to 1 m yr⁻¹; see Bea et al. [9] for details).

Locally, a decrease of solution density and specific ion concentrations can be observed in response to meltwater ingress during the period of glacial retreat (Figures 4(c) and 4(d)). Isolated lenses of fresh water reach depths of up to 350 m within confined aquifers that subcrop in the model domain (Figure 4(c)). Focusing on the upper 1,000 m of the basin and using chloride concentrations as an indicator for solution density, the ingress of glacial meltwater during the temperate glacial retreat period is observed (exemplified by comparing concentration distributions at = 17,000 and 20,000 years, Figure 5). Meltwater ingress to depths of approximately 350 m occurs principally within the more permeable dolostone and sandstone units present in the basin (Figure 5).

Figure 5: Spatial Cl concentration distributions at different times for the base case simulation (i.e., at 0, 17,000, 20,000, and 32,000 years), corresponding to the initial condition, the end of Stage II, Stage III, and the final condition.

In addition, the simulation results indicate that the basin remains generally thermally stable with limited perturbations of groundwater temperature in the shallow region as the ice advances (Figures 4(e) and 4(f)).

5.2. O₂ Ingress and Attenuation

In addition to the processes of dedolomitization, evaporite dissolution, and ion exchange (see Bea et al. [8, 9] for more discussion), redox-buffering mineral phases are also expected to contribute to water-rock interactions. In particular, the oxidative dissolution of reduced mineral phases and the oxidation of SOM are expected to significantly attenuate the depth of O₂ penetration in glacial recharge as is evident in the base case simulation during Stage III, the temperate glacier retreat (see Figure 3(a)).

The base case simulation demonstrates that the distribution and attenuation of O₂ is controlled by the presence and abundance of the O₂-consuming phases, hydrogeological processes, and the bed-boundary conditions assumed during the advance and retreat of the ice sheet. During quasi-steady-state interglacial conditions (e.g., = 0 years), only very limited ingress of O₂ into the sedimentary basin is seen, mostly focused on the dolostone units Dol2 and Dol3 (Figure 6). During the period of glacial advance and glacial maximum (i.e., Stages I and II), dissolved O₂ concentrations in the sedimentary basin are observed to decline relative to interglacial conditions. This behavior is due to the lack of recharge below the cold-based ice sheet (see results at 17,000 years in Figure 6). Oxygen present in the sedimentary units is consumed under these conditions by reactions with redox-buffering mineral phases and SOM (Figures 6 and 7). This trend is most clearly observed in the dolostone unit Dol2, where dissolved O₂ concentrations at 25 m depth abruptly decrease below the numerical detection limit as the ice sheet overruns this unit at = 4,000 years (Figure 7(c)). Oxygen does not reappear in this unit until the onset of glacial retreat and the establishment of warm-based glacial conditions. Similar behavior can be observed in other sedimentary units during the stage of glacial advance (Figure 7).

Figure 6: Spatial O₂(aq) concentration distributions at different times for the base case simulation (i.e., at 0, 17,000, 20,000, and 32,000 years), corresponding to the initial condition, the end of Stage II, Stage III, and the final condition. O₂ concentrations in mol l⁻¹.

Figure 7: Temporal evolution of dissolved O₂ and Cl concentrations at the different observation point locations for the base case simulation. (a), (b), (c), (d), (e), (f), and (g) for Dol3, Sand4, Dol2, Sand3, Dol1, Sand2, and Sand1, respectively. O₂ concentrations in mol l⁻¹. If O₂ profiles are not visible for a specific depth and unit, concentrations remain below 1 × 10⁻⁹ mol l⁻¹ for duration of simulation.

The maximum O₂ ingress takes place during glacial meltwater production and infiltration (i.e., Stage III) with simulation results showing the highest O₂-concentrations in the more permeable units (i.e., Sand1, Sand2, Sand3, Sand4, Dol2, and Dol3, see Figure 6, = 20,000 years). Oxygen ingress is also clearly depicted in the temporal concentration evolution at different observation locations (Figure 7). For instance, in Sand3 at 25 m depth, O₂ concentrations increase by almost five orders of magnitude during the retreat and melting stage; however, O₂ increase is negligible at greater depths (Figure 7(d)). Maximum O₂ concentrations at 25 m depth are already close to one order of magnitude lower than the O₂ content in the ingressing meltwater (i.e., 1.3 × 10⁻³ mol l⁻¹), demonstrating the important role of oxidative dissolution reactions.

The rapid O₂ increases seen at the onset of Stage III, mainly at 25 m depth in several of the more permeable units (see results at 17,500 years in Dol3, Sand3, Dol1, Sand2, and Sand1 in Figures 7(a), 7(d), 7(e), 7(f), and 7(g), resp.), can be attributed to rapid O₂ ingress as a consequence of the instantaneous change in boundary conditions to a warm-based ice sheet with an imposed piezometric head of 95% of the ice sheet thickness. The subsequent decline and increase in O₂ concentrations (Figure 7) results from the dynamic interactions between declining hydraulic heads as the ice sheet retreats and hydromechanical unloading, affecting hydraulic head distributions and groundwater velocities within the groundwater system.

Later during Stage III, O₂-concentrations diminish substantially in Dol3, Sand4, Dol2, and Sand3 and towards the end of this stage also in Sand2 and Sand1 as the ice sheet retreats. This decline is due to slowing recharge related to the decrease of the ice sheet thickness and lower hydraulic heads, limiting the driving force for meltwater ingress, while O₂-consuming reactions remain active at similar rates.

The most significant O₂ ingress is seen in the Dol2 and Dol3 units, even though these units do not possess the highest hydraulic conductivities. In these units, the relative lack of O₂ attenuation is due to the absence of the redox-buffering mineral phases and the low abundance of SOM (see Table 2). On the other hand, the more permeable sandstones show stronger attenuation potential, which can be attributed to the presence of pyrite (Table 2).

A comparison of the results for the depth of freshwater ingress, evidenced by the decline of chloride concentrations and the depth of O₂-ingress (compare Figures 5 and 6; see Figure 7), illustrates the enhanced attenuation of O₂. This outcome reaffirms the important role of oxidative dissolution reactions and SOM oxidation, inhibiting O₂ ingress in addition to hydraulic controls imposed by the architecture of the site-specific hydrostratigraphy and the presence of dense brines within the sedimentary basin groundwater system.

At the end of the simulation (representing interglacial recharge conditions), O₂ concentrations in the sedimentary units return to initial O₂ contents (Stage IV, see results at 32,000 years in Figure 6). This behavior is again most clearly observed at 25 m depth in the Dol2 unit (Figure 7(c)).

To obtain these results, reactive surface areas and the effective SOM oxidation rate were set to achieve a maximum O₂ ingress, at concentrations of 1 × 10⁻⁹ mol l⁻¹, to depths of 200 m during the period of deglaciation and all other times of the simulation. This maximum depth of ingress was achieved using surface areas of 5 × 10⁻³ m² g⁻¹ and an effective rate coefficient for SOM of 5 × 10⁻¹² mol dm⁻³ SOM s⁻¹ (Table 3). These reactive surface areas correspond to 25, 15, and 13 m² dm⁻³ mineral for pyrite, biotite, and chlorite, respectively.

The resulting reactive surface areas can be set in context with specific mineral surface areas reported in the literature. For instance, the specific surface area of framboidal pyrite is approximately 2 m² g⁻¹ for mean diameters ranging from 60 to 150 μm, with values as high as 4 m² g⁻¹ for diameters between 2 and 5 μm [40]. Similarly, the specific surface area for chlorite ranges between 3.4 and 4.4 m² g⁻¹ for grain sizes from 125 to 300 μm and 75 to 125 μm, respectively [41]. Beckingham et al. [42] report specific surface areas for biotite (0.5–4.7 m² g⁻¹, grain size range from 10 to 420 μm), chlorite (2.8–7.6 m² g⁻¹, 10–250 μm), and pyrite (0.03–1.1 m² g⁻¹, 10–250 μm). Although surface area is strongly dependent on mineral morphology and grain size, it can be concluded that the reactive surface areas considered for the base case are typically around three orders of magnitude lower than literature values for specific surface areas, implying that the base case mineral reactivities to limit O₂-ingress to 200 m are highly conservative, even taking into account the often large discrepancies between laboratory-measured specific surface areas and reactive surface areas in the field. It should be noted that these results were obtained using an equivalent porous medium approach; however, a large degree of conservatism remains, even if the majority of flow occurs within interconnected fracture networks or partially dolomitized horizons.

Using these parameters and the mineral volume fractions listed in Table 2, the maximum reaction rates in the various units can be calculated for the individual minerals and SOM (Table 6). Table 6 shows that the effective reaction rates for the base case range between 5 × 10⁻¹⁵ and 6 × 10⁻¹³ mol dm⁻³ bulk s⁻¹ in all units with the exception of the shales, which show higher reaction rates due to more abundant pyrite and chlorite. However, the permeability of the shale units is so low that the reactivity of redox-buffering mineral phases is irrelevant with respect to limiting O₂ ingress.

Table 6: Maximum oxidative dissolution rates and SOM oxidation rates for base case calculated based on parameters provided in Tables 2 and 3 and pO₂ = 1 atm, rates in units of mol dm⁻³ bulk s⁻¹.

Mineral
Cr
Crw
Sand1
Lim1
Sand2
Sh1
Dol1
Ev
Sh2
Sand3
Dol2
Sand4
Sh3
Dol3

Pyrite

Biotite

Chlorite

Organic matter

In comparison to literature data for oxidative mineral dissolution rates and SOM oxidation, effective rates on the order of 10⁻¹⁵ to 10-13 mol dm⁻³ bulk s⁻¹, as used for the base case, are very low. For example, utilizing rate data, mineral abundances, and effective surface areas reported by Beckingham et al. [42], effective dissolution rates for pyrite, biotite, and chlorite in a volcanic sandstone are determined as 2.3 × 10⁻⁹, 2.1 × 10⁻¹⁰, and 2.7 × 10⁻¹⁰ mol dm⁻³ bulk s⁻¹ for pyrite, biotite, and chlorite, respectively. Mineral volume fractions associated with these calculations are 0.4, 2.6, and 1.2 dm³ mineral dm⁻³ bulk for pyrite, biotite, and chlorite, respectively, within a factor of 5 of the volumetric fractions of the present study. Similarly, reaction rates for SOM reported in the literature for Pleistocene-Holocene sediments at 40 m depth [43] range from 7 × 10⁻⁷ to 4 × 10⁻⁹ mol dm⁻³ SOM s⁻¹ for fine and coarse fractions, respectively. These rates correspond to 7 × 10⁻¹⁰ to 4 × 10⁻¹² mol dm⁻³ bulk s⁻¹ for a low SOM volumetric fraction of 0.1% assumed here for the sandstone and dolostone units (Table 2). It can be concluded that typical reaction rates reported in the literature are at least two orders of magnitude faster than values used in the base case (Table 6).

Overall, these results indicate that the reaction parameters used for the base case are highly conservative (i.e., low). Even with this high level of conservatism, O₂ ingress is substantially inhibited by the presence of reduced mineral phases and SOM, even under conditions of glacial retreat and substantial O₂-rich meltwater production.

5.3. Sensitivity Analysis Results

To further investigate the role of O₂-consuming reactions and ingress during a glaciation cycle, a sensitivity analysis was conducted, including the following: (a) Scenario I, increasing the conservative base case reactions rates by one order of magnitude; (b) Scenario II, decreasing the O₂-concentration in the glacial meltwaters to less conservative atmospheric values; and (c) Scenario III, combining both Scenarios I and II.

The results demonstrate that an increase in effective reaction rates by one order of magnitude further limits O₂ ingress, even during periods of glacial meltwater ingress (Figures 8 and 9). With the increased rates, the maximum O₂ concentrations at 25 m depth are almost two orders of magnitude lower than those in the meltwater. The maximum O₂-ingress depth is now restricted to less than 100 m, with O₂ concentrations at 75 m slightly above the numerical detection limit only in the Dol2 and Sand3 units (see Figures 9(c) and 9(d)).

Figure 8: Spatial O₂(aq) concentration distributions at different times for Scenario 1 with a 10-fold increase of O₂-consumption rates (i.e., at 0, 17,000, 20,000, and 32,000 years), corresponding to the initial condition, the end of Stage II, Stage III, and the final condition. O₂ concentrations in mol l⁻¹.

Figure 9: Temporal evolution of dissolved O₂ and Cl concentrations at the different observation point locations for Scenario 1 with a 10-fold increase of O₂-consumption rates. (a), (b), (c), (d), (e), and (f) for Dol3, Sand4, Dol2, Sand3, Sand2, and Sand1, respectively. O₂ concentrations in mol l⁻¹. If O₂ profiles are not visible for a specific depth and unit, concentrations remain below 1 × 10⁻⁹ mol l⁻¹ for duration of simulation.

Comparison of the base case with the three sensitivity scenarios shows how a decline in O₂ in recharge waters and/or higher oxidative dissolution rates limit O₂ ingress at a depth of 75 m (Figure 10). For all three sensitivity scenarios, O₂ concentrations are restricted to the shallowest part of the aquifers, with lower penetration depths than the base case and with none of these scenarios showing O₂ concentrations above the numerical detection limit (1 × 10⁻⁹ mol l⁻¹, or 3.2 × 10-5 mg l⁻¹) at 150 m (results not shown).

Figure 10: Temporal evolution of dissolved O₂ concentrations at a depth of 75 m comparing results from Scenarios I–III to the results of the base case. (a), (b), (c), (d), (e), and (f) for Dol3, Sand4, Dol2, Sand3, Sand2, and Sand1, respectively. O₂ concentrations in mol l⁻¹. If O₂ profiles are not visible for a specific unit, concentrations remain below 1 × 10⁻⁹ mol l⁻¹ for duration of simulation.

6. Summary and Concluding Remarks

Oxygen migration and consumption in a large-scale hypothetical intracratonic sedimentary basin during a 32,500-year glaciation-deglaciation cycle were investigated using reactive transport modeling. The sensitivity of O₂-ingress to O₂ concentration in recharge water and rates of O₂-consuming reactions was evaluated.

Although illustrative in nature, the numerical model results were intended to capture the main hydrogeological features of the Paleozoic sedimentary basins in North America. To the degree possible, model parameters were based on published data or relationships.

Analysis of a base case reactive transport simulation in this hypothetical sedimentary basin suggests that the evolution of hydraulic head distributions due to mechanical loading during the period of glacial advance and meltwater production during the period of glacial retreat may result in substantial hydraulic gradients. However, the hydrostratigraphy of the interbedded low and high hydraulic conductivity sedimentary rock formations significantly limit groundwater flow velocities and mixing between aquifers and aquitards. The presence of dense brines at depth exerts an important control on fresh water ingress through the maintenance of density gradients, despite the transient hydraulic and mechanical surface boundary conditions created by an advancing and retreating ice sheet.

The key objective of the study was to assess the role of mineral oxidative dissolution reactions and SOM oxidation in inhibiting O₂-ingress into sedimentary basins. While the simulations suggest that freshwater ingress can occur to a maximum depth of up to 350 m during the deglaciation period, it was shown that even for highly conservative assumptions regarding O₂ consumption rates, O₂ ingress was limited to depths of less than 200 m at concentrations not exceeding 1 × 10⁻⁹ mol l⁻¹. Rates used to allow O₂ to reach depths of 200 m for the mineral phases and SOM were typically 3–5 orders of magnitude lower than rates reported in the literature, implying a high level of conservatism. For slightly higher, yet still conservative rates and/or lower O₂ concentrations in glacial meltwater, O₂ ingress was even lower, essentially restricted to the top 100 m of the sedimentary basin. The simulation results suggest that detailed characterization of the present-day abundance and spatial distributions of reduced mineral phases and SOM would be beneficial to assess O₂ attenuation capacity of sedimentary basins. On the other hand, it appears that the difficult task of characterizing field-scale effective O₂ consumption rates might not be necessary, considering that limited O₂-ingress was achieved based on conservative assumptions.

However, the highly complex nature of the simulated system must be taken into consideration and the results must be viewed with care. Although the simulations were constrained by field observations to the degree possible, assumptions had to be made regarding several parameters and processes, including the abundance and distribution of reactive mineral phases and SOM, recharge patterns during periods of glaciation and deglaciation, and the hydraulic and mechanical properties of the various sedimentary units. In addition, simplifications were made regarding dimensionality (restricted to two spatial dimensions), the description of flow and transport processes (an Equivalent Porous Medium approach was used, neglecting the presence of fracture networks), mineralogical composition of the rocks, and the form of the rate expressions describing oxidative dissolution reactions and SOM oxidation. These assumptions and simplifications lead to substantial uncertainties associated with these simulations and limit their predictive capability.

Nevertheless, overall, our results suggest a high degree of geochemical stability of sedimentary basins and a significant capacity to attenuate O₂ in recharge water and to maintain reducing conditions at depth over long time periods including during glaciation-deglaciation events.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

Funding for this research was provided by a NWMO (Nuclear Waste Management Organization, Canada) grant held by K. Ulrich Mayer and Kerry T. B. MacQuarrie. The authors would like to acknowledge WestGrid and Compute Canada for providing computing hardware, software, and technical support. Also they are grateful for financial support from CONICET (National Council of Scientific and Technical Research), Argentina.

References

M. Mazurek, Long-term used fuel waste management-Geoscientific review of the sedimentary sequence in southern Ontario. Background Paper 612 prepared by Rock-Water Interactions, Institute of Geological Sciences, University of Bern for the Nuclear Waste Management Organization, 2004.

S. M. Spiessl, K. T. B. MacQuarrie, and K. U. Mayer, “Identification of key parameters controlling dissolved oxygen migration and attenuation in fractured crystalline rocks,” Journal of Contaminant Hydrology, vol. 95, no. 3-4, pp. 141–153, 2008.

M. Sidborn and I. Neretnieks, “Long-term oxygen depletion from infiltrating groundwaters: Model development and application to intra-glaciation and glaciation conditions,” Journal of Contaminant Hydrology, vol. 100, no. 1-2, pp. 72–89, 2008.

J. C. McIntosh and L. M. Walter, “Paleowaters in Silurian-Devonian carbonate aquifers: Geochemical evolution of groundwater in the Great Lakes region since the Late Pleistocene,” Geochimica et Cosmochimica Acta, vol. 70, no. 10, pp. 2454–2479, 2006.

J. C. McIntosh, L. M. Walter, and A. M. Martini, “Pleistocene recharge to midcontinent basins: Effects on salinity structure and microbial gas generation,” Geochimica et Cosmochimica Acta, vol. 66, no. 10, pp. 1681–1700, 2002.

J. C. McIntosh and L. M. Walter, “Volumetrically significant recharge of Pleistocene glacial meltwaters into epicratonic basins: constraints imposed by solute mass balances,” Chemical Geology, vol. 222, no. 3-4, pp. 292–309, 2005.

J. C. McIntosh, G. Garven, and J. S. Hanor, “Impacts of Pleistocene glaciation on large-scale groundwater flow and salinity in the Michigan Basin,” Geofluids, vol. 11, no. 1, pp. 18–33, 2011.

S. Bea, K. U. Mayer, and K. T. B. MacQuarrie, Modelling Reactive Transport in Sedimentary Rock Environments-Phase II, Reactive Transport Modeling in a Sedimentary Basin Affected by a Glaciation/Deglaciation Event, Nuclear Waste Management Organization, 2011.

S. A. Bea, U. K. Mayer, and K. T. B. MacQuarrie, “Reactive transport and thermo-hydro-mechanical coupling in deep sedimentary basins affected by glaciation cycles: model development, verification, and illustrative example,” Geofluids, vol. 16, no. 2, pp. 279–300, 2016.

A. M. Martini, L. M. Walter, J. M. Budat, T. C. W. Ku, C. J. Kaiser, and M. Schoell, “Genetic and temporal relations between formation waters and biogenic methane: Upper Devonian Antrim shale, Michigan Basin, USA,” Geochimica et Cosmochimica Acta, vol. 62, no. 10, pp. 1699–1720, 1998.

J. C. McIntosh, L. M. Walter, and A. M. Martini, “Extensive microbial modification of formation water geochemistry: Case study from a Midcontinent sedimentary basin, United States,” Bulletin of the Geological Society of America, vol. 116, no. 5-6, pp. 743–759, 2004.

T. R. Weaver, S. K. Frape, and J. A. Cherry, “Recent cross-formational fluid flow and mixing in the shallow Michigan Basin,” Bulletin of the Geological Society of America, vol. 107, no. 6, pp. 697–707, 1995.

J. J. Kolak, D. T. Long, J. M. Matty, G. J. Larson, D. F. Sibley, and T. B. Councell, “Ground-water, large-lake interactions in Saginaw Bay, Lake Huron: A geochemical and isotopic approach,” Bulletin of the Geological Society of America, vol. 111, no. 2, pp. 177–188, 1999.

M. M. Husain, J. A. Cherry, and S. K. Frape, “The persistence of a large stagnation zone in a developed regional aquifer, southwestern Ontario,” Canadian Geotechnical Journal, vol. 41, no. 5, pp. 943–958, 2004.

L. Ma, M. C. Castro, and C. M. Hall, “A late Pleistocene-Holocene noble gas paleotemperature record in southern Michigan,” Geophysical Research Letters, vol. 31, no. 23, pp. 1–4, 2004.

S. Klump, T. Grundl, R. Purtschert, and R. Kipfer, “Groundwater and climate dynamics derived from noble gas, 14C, and stable isotope data,” Geology, vol. 36, no. 5, pp. 395–398, 2008.

T. Wolery and S. Daveler, EQ6, A computer program for reaction path modeling of aqueous geochemical system: Theoretical manual, user's guide, and related documentation (version 7.0). UCLR-MA-110662 PT IV, Lawrence Livermore Natl. Lab. Livermore, California, 1992.

C. Appelo and D. Postma, Geochemistry, Groundwater and Pollution, A. A. Balkema, Ed., 1993.

M. Hobbs, S. Frape, O. Shouakar-Stash, and L. Kennell, Phase I Regional Hydrogeochemistry, Southern Ontario. Ontario Power Generation, Department of Earth and Environmental Science, University of Waterloo, OPG 00216-REP-01300-00006-R00, 2008.

J. Guimerà, L. Duro, and A. Delos, Changes in groundwater composition as a consequence of deglaciation: implications for performance assessment. Report no. R-06-105. Svensk Kärnbränslehantering AB, 2006.

K. U. Mayer, E. O. Frind, and D. W. Blowes, “Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions,” Water Resources Research, vol. 38, no. 9, pp. 13-1–13-21, 2002.

S. A. Bea, S. A. Wilson, K. U. Mayer, G. M. Dipple, I. M. Power, and P. Gamazo, “Reactive transport modeling of natural carbon sequestration in ultramafic mine tailings,” Vadose Zone Journal, vol. 11, no. 2, 2012.

D. Su, K. U. Mayer, and K. T. B. MacQuarrie, “Parallelization of MIN3P-THCm: A high performance computational framework for subsurface flow and reactive transport simulation,” Environmental Modeling and Software, vol. 95, pp. 271–289, 2017.

V. F. Bense and M. A. Person, “Transient hydrodynamics within intercratonic sedimentary basins during glacial cycles,” Journal of Geophysical Research: Earth Surface, vol. 113, no. 4, 2008.

S. E. Grasby and Z. Chen, “Subglacial recharge into the Western Canada Sedimentary Basin—impact of Pleistocene glaciation on basin hydrodynamics,” Bulletin of the Geological Society of America, vol. 117, no. 3-4, pp. 500–514, 2005.

A. M. Provost, C. I. Voss, and C. E. Neuzil, “Glaciation and regional groundwater flow in the Fennoscandian shield,” Geofluids, vol. 12, no. 1, pp. 79–96, 2012.

C. E. Neuzil, “Hydromechanical coupling in geologic processes,” Hydrogeology Journal, vol. 11, no. 1, pp. 41–83, 2003.

T. H. Henderson, K. U. Mayer, B. L. Parker, and T. A. Al, “Three-dimensional density-dependent flow and multicomponent reactive transport modeling of chlorinated solvent oxidation by potassium permanganate,” Journal of Contaminant Hydrology, vol. 106, no. 3-4, pp. 195–211, 2009.

C. E. Harvie, N. Møller, and J. H. Weare, “The prediction of mineral solubilities in natural waters: The Na-K-Mg-Ca-H-Cl-SO₄-OH-HCO₃-CO₃-CO₂-H₂O system to high ionic strengths at 25°C,” Geochimica et Cosmochimica Acta, vol. 48, no. 4, pp. 723–751, 1984.

K. S. Pitzer, “Thermodynamics of electrolytes. I. theoretical basis and general equations,” The Journal of Physical Chemistry C, vol. 77, no. 2, pp. 268–277, 1973.

C. Monnin, “Density calculation and concentration scale conversions for natural waters,” Computers & Geosciences, vol. 20, no. 10, pp. 1435–1445, 1994.

C. Voss and A. Provost, “SUTRA-A Model for saturated-unsaturated variable-density ground-water flow with solute or energy transport (Version 2.1), 2008”.

D. B. Bahr, E. W. H. Hutton, J. P. M. Syvitski, and L. F. Pratson, “Exponential approximations to compacted sediment porosity profiles,” Computers & Geosciences, vol. 27, no. 6, pp. 691–700, 2001.

S. T. Brennan and T. K. Lowenstein, “The major-ion composition of Silurian seawater,” Geochimica et Cosmochimica Acta, vol. 66, no. 15, pp. 2683–2700, 2002.

J. S. Hanor and J. C. McIntosh, “Are secular variations in seawater chemistry reflected in the compositions of basinal brines?” Journal of Geochemical Exploration, vol. 89, no. 1-3, pp. 153–156, 2006.

USDOE (United States Department of Energy): In-drift precipitates/salts model. Argonne National Laboratory Report EBS-MD-000045 REV 03. Chicago, USA, 2007.

S. A. Bea, J. Carrera, C. Ayora, and F. Batlle, “Modeling of concentrated aqueous solutions: Efficient implementation of Pitzer equations in geochemical and reactive transport models,” Computers & Geosciences, vol. 36, no. 4, pp. 526–538, 2010.

M. A. Williamson and J. D. Rimstidt, “The kinetics and electrochemical rate-determining step of aqueous pyrite oxidation,” Geochimica et Cosmochimica Acta, vol. 58, no. 24, pp. 5443–5454, 1994.

J. Johnson and J. L. Fastook, “Northern Hemisphere glaciation and its sensitivity to basal melt water,” Quaternary International, vol. 95, no. 96, pp. 65–74, 2002.

C. E. Pugh, L. R. Hossner, and J. B. Dixon, “Pyrite and marcasite surface area as influenced by morphology and particle diameter.,” Soil Science Society of America Journal, vol. 45, no. 5, pp. 979–982, 1981.

M. Malmström, S. Banwart, J. Lewenhagen, L. Duro, and J. Bruno, “The dissolution of biotite and chlorite at 25°C in the near-neutral pH region,” Journal of Contaminant Hydrology, vol. 21, no. 1-4, pp. 201–213, 1996.

L. E. Beckingham, E. H. Mitnick, C. I. Steefel et al., “Evaluation of mineral reactive surface area estimates for prediction of reactivity of a multi-mineral sediment,” Geochimica et Cosmochimica Acta, vol. 188, pp. 310–329, 2016.

N. Hartog, J. Griffioen, and C. H. Van der Weijden, “Distribution and reactivity of O2-reducing components in sediments from a layered aquifer,” Environmental Science & Technology, vol. 36, no. 11, pp. 2338–2344, 2002.