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

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-MIN3PTHCm 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 and 6 × 10mol dm bulk s were sufficient to restrict the depth of O2 ingress to less than 200m.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.


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 2 ) 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 2 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][5][6][7][8][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 14 C ages correspond to periods of ice cover, and the low  18 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, 14 C ages ranging from modern to >50 ka, and oxygen isotope values ranging from −18 to −13‰ [4,[12][13][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 2 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 2 -concentrations [2,20].It is therefore of interest to evaluate the probability for enhanced dissolved O 2 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 2 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 2 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 2 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 2 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 2 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 2 , 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][22][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 2 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 2 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 2 consumption rates in sedimentary basins, the model is used to identify the effective reaction rates required to restrict O 2 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 2 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 2 ingress towards effective rates and O 2 content in recharge water is investigated.

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][5][6][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 4 ], 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 2 and NaCltype), 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 regionalscale 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][25][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 lowpermeability units like shales (e.g., [8,24,27]).These effects on the flow regime may also affect the ingress and fate of dissolved O 2 into the sedimentary units.
Under present-day conditions, advective solute migration in shallow permeable aquifer units, including transport of dissolved O 2 contained in the surficial recharge water, is limited to shallow depths [8,9].In addition, dissolved O 2 could potentially be consumed by reactions with reduced mineral phases (e.g., pyrite, chlorite, or biotite) and SOM, further contributing to O 2 -depleted conditions near the ground surface.The depth to which O 2 can migrate will depend on the interplay between altered fluid circulation patterns, the O 2 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.

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 (1) 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; (2) fluid density calculations based on fluid volumetric predictions derived from the Pitzer equations according to Monnin [31]; (3) one-dimensional hydromechanical coupling [27] due to ice sheet advance and retreat; and (4) 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].

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]).

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 4 , K, Total Inorganic Carbon (TIC), Na, Br, O 2 , and H 2 O.For the sake of simplicity, the components Al, Fe(II), Fe(III), and H 4 SiO 4 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]   stoichiometrically as CH 2 O).These redox-buffering solids were chosen based on their typical abundance and level of reactivity.Pyrite reacts readily with O 2 (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 O 2 .The rate expression for the oxidative dissolution of pyrite is derived from the work of Weaver et al. [12], assuming that only the O 2 pathway is active, with a square root dependence on O 2 .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  pyr [mol dm −3 bulk s −1 ] is the oxidative dissolution rate of pyrite, [O 2 ] is the dissolved oxygen concentration [mol l −1 ],  pyr is the kinetic rate constant evaluated at pH = 7 [mol 0.5 l 0.5 m −2 mineral s −1 ],  pyr is the mineral volumetric fraction of pyrite [dm 3 mineral dm −3 bulk],  pyr is the density of pyrite [g dm −3 mineral], and 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 [dm 3 mineral dm −3 bulk],   is the relevant mineral density [g dm −3 mineral], and 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 (1) (1) (1) Equilibrium constant taken from Wolery and Daveler [17]. (2)Selectivity coefficients taken Appelo and Postma [18]. (3)Not controlled by equilibrium [2].
distribution of organic matter in the different sedimentary units: ) , where and  SOM is the volumetric fraction of organic matter present in the solid phase [dm 3 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 O 2 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 2 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 O 2 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 Table 2: Initial mineral contents (%) considered in the reactive transport simulations (for reference see Bea et al., [9]).Based on rock composition for Queenston, Georgian Bay, and Blue Mountain Formations reported by Mazurek [1].Table 3: Kinetic parameters and solid properties for the redox-buffering minerals used in base case simulation. (1)Redox-buffering mineral

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]): (1) Stage I corresponds to a phase of ice sheet accumulation (i.e., loading, for a period of 12,500 years), (2) Stage II represents a constant ice sheet thickness for a period of 5,000 years, (3) Stage III corresponds to the phase of ice sheet melting and retreat (i.e., unloading, for a period of 5,000 years), and (4) Stage IV is representative of conditions with no ice present (i.e., return to present-day surface boundary conditions).
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 O 2 ingress occurs beneath the ice sheet because cold-based conditions were assumed.During glacial retreat, elevated O 2 concentrations were assumed in the meltwater beneath the ice sheet.

Evaluation of O 2 Ingress and Attenuation
To investigate the impact of a glaciation cycle and uncertain geochemical parameters on O 2 ingress and consumption, a   ,  eff,SOM increased 10-fold 0.21 (1) Partial pressure of O 2 (g) in atm corresponds to 60 mg l −1 for 1 atm and 12.6 mg l −1 for 0.21 atm.
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 2 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 2 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 2 consumption rates and O 2 concentrations in source glacial meltwater on depth of O 2 penetration.

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 2 attenuation and redox stability of the basin.

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

O 2 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 waterrock interactions.In particular, the oxidative dissolution of reduced mineral phases and the oxidation of SOM are expected to significantly attenuate the depth of O 2 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 2 is controlled by the presence and abundance of the O 2 -consuming phases, hydrogeological processes, and the bed-boundary conditions assumed during the advance and retreat of the ice sheet.During quasisteady-state interglacial conditions (e.g.,  = 0 years), only very limited ingress of O 2 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 2 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 redoxbuffering mineral phases and SOM (Figures 6 and 7).This trend is most clearly observed in the dolostone unit Dol2, where dissolved O 2 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 O 2 ingress takes place during glacial meltwater production and infiltration (i.e., Stage III) with simulation results showing the highest O 2 -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 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 2 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 2 -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 2 -consuming reactions remain active at similar rates.
The most significant O 2 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 2 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 2 -ingress (compare Figures 5 and 6; see Figure 7), illustrates the enhanced attenuation of O 2 .This outcome reaffirms the important role of oxidative dissolution reactions and SOM oxidation, inhibiting O 2 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 2 concentrations in the sedimentary units return to initial O 2 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 2 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 m 2 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 m 2 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 m 2 g −1 for mean diameters ranging from 60 to 150 m, with values as high as 4 m 2 g −1 for diameters between 2 and 5 m [40].Similarly, the specific surface area for chlorite ranges between 3.4 and 4.4 m 2 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 m 2 g −1 , grain size range from 10 to 420 m), chlorite (2.8-7.6 m 2 g −1 , 10-250 m), and pyrite (0.03-1.1 m 2 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 O 2 -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 O 2 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 dm 3 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, O 2 ingress is substantially inhibited by the presence of reduced mineral phases and SOM, even under conditions of glacial retreat and substantial O 2 -rich meltwater production.

Sensitivity Analysis Results.
To further investigate the role of O 2 -consuming reactions and ingress during a glaciation cycle, a sensitivity analysis was conducted, including the following: (a) Scenario I, increasing the conservative base 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 The results demonstrate that an increase in effective reaction rates by one order of magnitude further limits O 2 ingress, even during periods of glacial meltwater ingress (Figures 8 and 9).With the increased rates, the maximum O 2 concentrations at 25 m depth are almost two orders of magnitude lower than those in the meltwater.The maximum O 2 -ingress depth is now restricted to less than 100 m, with O 2 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 O 2 in recharge waters and/or higher oxidative dissolution rates limit O 2 ingress at a depth of 75 m (Figure 10).For all three sensitivity scenarios, O 2 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 2 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).

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 2 -ingress to O 2 concentration in recharge water and rates of O 2 -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 2 -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 2 consumption rates, O 2 ingress was 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 2 attenuation capacity of sedimentary basins.On the other hand, it appears that the difficult task of characterizing field-scale effective O 2 consumption rates might not be necessary, considering that limited O 2 -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 2 in recharge water and to maintain reducing conditions at depth over long time periods including during glaciation-deglaciation events.

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

1 N
o n r e a c t i v em i n e r a l
Evolution of maximum ice sheet thickness and extension topography-driven flow ice sheet (cold-based) Boundary conditions, Stages I and II Impervious boundary topography-driven flow ice sheet (warm-based) retreats water table Boundary conditions, Stage III topography-driven flow

Figure 4 :
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).

Figure 5 :Figure 6 : 1 )Figure 7 :
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.

2 = 1 Figure 8 :
Figure 8: Spatial O 2 (aq) concentration distributions at different times for Scenario 1 with a 10-fold increase of O 2 -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 2 concentrations in mol l −1 .

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

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

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