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 (O2). 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, O2-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 O2, 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−15 and 6 × 10−13 mol dm−3 bulk s−1 were sufficient to restrict the depth of O2 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 O2 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 (O2) 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 O2 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., [49]). 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 14C ages correspond to periods of ice cover, and the low δ18O 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, 14C ages ranging from modern to >50 ka, and oxygen isotope values ranging from −18 to −13 [4, 1214]. 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 O2 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 O2-concentrations [2, 20]. It is therefore of interest to evaluate the probability for enhanced dissolved O2 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 O2 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 O2 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 O2 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 O2 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 O2 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 O2, 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, 2123] 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 O2 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 O2 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 O2 consumption rates in sedimentary basins, the model is used to identify the effective reaction rates required to restrict O2 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 O2 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 O2 ingress towards effective rates and O2 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, 47]). 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 [CaSO4], 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 (CaCl2 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, 2426]. 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 O2 into the sedimentary units.

Under present-day conditions, advective solute migration in shallow permeable aquifer units, including transport of dissolved O2 contained in the surficial recharge water, is limited to shallow depths [8, 9]. In addition, dissolved O2 could potentially be consumed by reactions with reduced mineral phases (e.g., pyrite, chlorite, or biotite) and SOM, further contributing to O2-depleted conditions near the ground surface. The depth to which O2 can migrate will depend on the interplay between altered fluid circulation patterns, the O2 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.

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, SO4, K, Total Inorganic Carbon (TIC), Na, Br, O2, and H2O. For the sake of simplicity, the components Al, Fe(II), Fe(III), and H4SiO4 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.

Mineral phases considered include halite [NaCl], anhydrite [CaSO4], calcite [CaCO3], and dolomite [CaMg(CO3)2], together with the redox-buffering solids pyrite [FeS2], biotite [K(Mg2Fe)(Si3Al)O10(OH)2], chlorite [Mg2Fe3Al2Si3O10(OH)8], and SOM (represented stoichiometrically as CH2O). These redox-buffering solids were chosen based on their typical abundance and level of reactivity. Pyrite reacts readily with O2 (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)).

The aforementioned redox-buffering mineral phases (i.e., pyrite, biotite, and chlorite) and SOM were modeled as irreversible kinetic reactions in the presence of O2. The rate expression for the oxidative dissolution of pyrite is derived from the work of Weaver et al. [12], assuming that only the O2 pathway is active, with a square root dependence on O2. 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 yieldwhere [mol dm−3 bulk s−1] is the oxidative dissolution rate of pyrite, is the dissolved oxygen concentration [mol l−1], is the kinetic rate constant evaluated at pH = 7 [mol0.5 l0.5 m−2 mineral s−1], is the mineral volumetric fraction of pyrite [dm3 mineral dm−3 bulk], is the density of pyrite [g dm−3 mineral], and is reactive surface area of pyrite [m2 g−1].

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−3 bulk s−1] is the oxidative dissolution rate of chlorite or biotite, is the corresponding kinetic rate constant [mol m−2 mineral s−1], is the mineral volumetric fraction of the phase considered [dm3 mineral dm−3 bulk], is the relevant mineral density [g dm−3 mineral], and is the reactive mineral surface area [m2 g−1].

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−3 bulk s−1] is the oxidation rate of SOM, is the effective rate coefficient [mol dm−3 SOM s−1], and is the volumetric fraction of organic matter present in the solid phase [dm3 SOM dm−3 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 O2 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 O2 ingress with concentrations above 1 × 10−9 mol l−1 (corresponding to 3.2 × 10−5 mg l−1) to a depth of 200 m. As indicated above, this depth was selected to conservatively define a tolerable depth of O2 ingress with regard to deep geologic repositories.

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.

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−2 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).

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 O2 ingress occurs beneath the ice sheet because cold-based conditions were assumed. During glacial retreat, elevated O2 concentrations were assumed in the meltwater beneath the ice sheet.

4. Evaluation of O2 Ingress and Attenuation

To investigate the impact of a glaciation cycle and uncertain geochemical parameters on O2 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 O2 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 O2 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 O2 consumption rates and O2 concentrations in source glacial meltwater on depth of O2 penetration.

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 O2 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.

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−1 to 1 m yr−1; 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).

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. O2 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 O2 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 O2 is controlled by the presence and abundance of the O2-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 O2 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 O2 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 O2 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).

The maximum O2 ingress takes place during glacial meltwater production and infiltration (i.e., Stage III) with simulation results showing the highest O2-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, O2 concentrations increase by almost five orders of magnitude during the retreat and melting stage; however, O2 increase is negligible at greater depths (Figure 7(d)). Maximum O2 concentrations at 25 m depth are already close to one order of magnitude lower than the O2 content in the ingressing meltwater (i.e., 1.3 × 10−3 mol l−1), demonstrating the important role of oxidative dissolution reactions.

The rapid O2 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 O2 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 O2 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, O2-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 O2-consuming reactions remain active at similar rates.

The most significant O2 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 O2 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 O2-ingress (compare Figures 5 and 6; see Figure 7), illustrates the enhanced attenuation of O2. This outcome reaffirms the important role of oxidative dissolution reactions and SOM oxidation, inhibiting O2 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), O2 concentrations in the sedimentary units return to initial O2 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 O2 ingress, at concentrations of 1 × 10−9 mol l−1, 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−3 m2 g−1 and an effective rate coefficient for SOM of 5 × 10−12 mol dm−3 SOM s−1 (Table 3). These reactive surface areas correspond to 25, 15, and 13 m2 dm−3 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 m2 g−1 for mean diameters ranging from 60 to 150 μm, with values as high as 4 m2 g−1 for diameters between 2 and 5 μm [40]. Similarly, the specific surface area for chlorite ranges between 3.4 and 4.4 m2 g−1 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 m2 g−1, grain size range from 10 to 420 μm), chlorite (2.8–7.6 m2 g−1, 10–250 μm), and pyrite (0.03–1.1 m2 g−1, 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 O2-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−15 and 6 × 10−13 mol dm−3 bulk s−1 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 O2 ingress.

In comparison to literature data for oxidative mineral dissolution rates and SOM oxidation, effective rates on the order of 10−15 to 10-13 mol dm−3 bulk s−1, 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−9, 2.1 × 10−10, and 2.7 × 10−10 mol dm−3 bulk s−1 for pyrite, biotite, and chlorite, respectively. Mineral volume fractions associated with these calculations are 0.4, 2.6, and 1.2 dm3 mineral dm−3 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−7 to 4 × 10−9 mol dm−3 SOM s−1 for fine and coarse fractions, respectively. These rates correspond to 7 × 10−10 to 4 × 10−12 mol dm−3 bulk s−1 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, O2 ingress is substantially inhibited by the presence of reduced mineral phases and SOM, even under conditions of glacial retreat and substantial O2-rich meltwater production.

5.3. Sensitivity Analysis Results

To further investigate the role of O2-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 O2-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 O2 ingress, even during periods of glacial meltwater ingress (Figures 8 and 9). With the increased rates, the maximum O2 concentrations at 25 m depth are almost two orders of magnitude lower than those in the meltwater. The maximum O2-ingress depth is now restricted to less than 100 m, with O2 concentrations at 75 m slightly above the numerical detection limit only in the Dol2 and Sand3 units (see Figures 9(c) and 9(d)).

Comparison of the base case with the three sensitivity scenarios shows how a decline in O2 in recharge waters and/or higher oxidative dissolution rates limit O2 ingress at a depth of 75 m (Figure 10). For all three sensitivity scenarios, O2 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 O2 concentrations above the numerical detection limit (1 × 10−9 mol l−1, or 3.2 × 10-5 mg l−1) at 150 m (results not shown).

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 O2-ingress to O2 concentration in recharge water and rates of O2-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 O2-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 O2 consumption rates, O2 ingress was limited to depths of less than 200 m at concentrations not exceeding 1 × 10−9 mol l−1. Rates used to allow O2 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 O2 concentrations in glacial meltwater, O2 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 O2 attenuation capacity of sedimentary basins. On the other hand, it appears that the difficult task of characterizing field-scale effective O2 consumption rates might not be necessary, considering that limited O2-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 O2 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.