Abstract
Geoscientific investigations for a proposed deep geologic repository at the Bruce Site, located on the eastern flank of the Michigan Basin, have identified unique and significant underpressured conditions. Along with the measurement of environmental tracer profiles (e.g., helium), this study aims to explore, through a series of numerical simulations, the nature of long-term phenomena responsible for the generation and preservation of formation underpressures. Three families of inverse numerical experiments for underpressure formation were examined by means of one-dimensional hydromechanically coupled models through the vertical hydrostratigraphic column: (i) uncertainty in glaciation scenarios; (ii) uncertainty in initial heads prior to glaciation; and (iii) uncertainty in the degree of hydraulic connectivity between the more permeable Guelph Formation at the Bruce Site and the applied glacial loading, for a total of 20 scenarios, assuming fully saturated conditions. Underpressured initial heads for the paleohydrogeologic simulations lead to lower calibrated vertical hydraulic conductivities. The robustness and resilience of the groundwater system to external perturbations are greater for the state where underpressured conditions predate the onset of glaciation and are better able to preserve the present day helium tracer profile in 260 Ma exhumation analyses.
1. Introduction
A deep geologic repository (DGR) for low and intermediate level radioactive wastes has been proposed at the Bruce Site in the Municipality of Kincardine in Ontario, Canada. The Bruce Site is located on the eastern flank of the intracratonic Paleozoic age Michigan Basin. A geologic cross-section through the basin indicating the location of the Bruce Site is shown in Figure 1. A diagram showing the site-specific sedimentary bedrock stratigraphy above the Precambrian basement that is comprised of 34 bedrock formations of Cambrian to Devonian age is illustrated in Figure 2.
Site characterization activities between 2006 and 2010 at the Bruce Site provided detailed information on the geologic, geochemical, hydrogeologic, and geomechanical properties of the sedimentary sequence in support of a geosynthesis and an on-going environmental assessment. In particular, hydraulic formation pressure data obtained during in situ hydraulic testing [1] and monitoring of three multilevel instrumented boreholes (DGR-series) through the sedimentary sequence indicated the occurrence of abnormal fluid pressures. Hydraulic heads within the Ordovician sequence were as much as 300 m below hydrostatic and the underlying permeable Cambrian sandstone was overpressured (Figure 3). Sedimentary basins, when at hydrological equilibrium, normally show a near-hydrostatic pressure distribution [2]. Under certain imposed conditions, transient abnormal pressures can develop in low permeability layers or other hydraulically isolated parts of groundwater systems [3].
Environmental and isotopic tracers are essential for characterizing low permeability settings in terms of fluid movement and residence times. These tracers represent natural analogues that provide a means to determine mass transport characteristics of solutes (e.g., radionuclides) over geologic time scales [6–10]. In various studies, helium isotopes have been used to quantify fluxes in low permeability media and to determine media properties [6, 11–15]. Helium in groundwater originates from the natural radioactive decay of uranium and thorium and from degassing of the crust and mantle [16, 17].
At the Bruce Site, the Ordovician shale and carbonate formations have preserved brine compositions with minimal evidence of transport since the Silurian, with this being attributed to permeability reduction resulting from secondary halite mineralization [9]. Clark et al. [9] determined helium isotope ratios with depth at the Bruce Site (see Figure 4) which show a distinct increase at the base of the Cobourg formation and further suggest that the Ordovician behaves as an aquiclude with limited transport capability over Paleozoic time. Xiang et al. [18] have measured effective diffusion coefficients at the Bruce Site on the order of 10−12 to 10−14 m2/s, indicative of extremely low solute transport. Numerical simulations by Al et al. [10] suggest that the vertical effective diffusion coefficients of the Ordovician carbonates are up to an order-of-magnitude lower at the formation scale than laboratory estimates. A comparison of Bruce Site horizontal hydraulic conductivities and effective diffusion coefficients with international data reveals that both key solute transport parameters are among the lowest recorded anywhere [4, 19, 20].
The Bruce Site investigations provide an opportunity, given the measurement of the unique and significant underpressured conditions and geochemical profiles, to examine the phenomena responsible for their occurrence and for estimating upscaled formation hydrogeologic properties. It further provides a compelling self-analogue that enables characterization at the Bruce Site of deep-seated groundwater system evolution as influenced by processes occurring on geologic time and space scales not otherwise achievable. A self-analogue refers to past geologic behaviour at a site which can be used as an analogue of the future behaviour expected to cause the transport of radionuclides at the site [21].
Numerous hypotheses regarding the phenomena responsible, solely or in combination, for the formation of hydraulic underpressures in sedimentary rock can be found in the literature. They include, among other osmosis, the presence of a nonwetting gas phase, hydromechanical effects due to exhumation, cyclic glacial loading, and crustal flexure [22–30]. Abnormal underpressures in geologic formations have been reported in many regions of the globe [2], including China [31, 32], Canada [33–35], and Switzerland [36]. Such abnormal pressures cannot exist without low permeability formations [3, 31, 36].
Low permeability sediments such as clays and shales may exhibit membrane behaviour resulting in the restriction of solute transport relative to the flow of water. This property gives rise to chemical osmosis: fluid flow in response to a chemical gradient [23]. In geologic media with membrane-like properties, osmotic behaviour can drive fluid flow from regions of low to regions of high solute concentrations, thereby raising fluid pressures where solute concentrations are high and decreasing fluid pressures where concentrations are low. The actual role of osmosis, however, with respect to influencing fluid pressure distributions has been controversial as the membrane properties of geologic media are not well understood [22]. Within the Pierre Shale in South Dakota, USA, in situ fluid pressure and solute concentration measurement over a nine-year period were found to be consistent with the generation of large (up to 20 MPa) osmotic pressure anomalies, requiring small shale porosity and large TDS contrasts [22]. At the Bruce Site, the underpressures observed in the DGR boreholes occur in both the Ordovician shale and carbonate sequences [4]. Carbonates do not function as membranes in the same manner as shales. Additionally, arguing against an osmotic explanation of the underpressures is the estimation of significant underpressures in the short-term (1–3 days) in situ straddle packer hydraulic tests in the DGR boreholes [37]; pressure changes as a result of osmosis would be expected to evolve at a much slower rate (hundreds to thousands of days) as demonstrated through field experiments by Neuzil [22].
Vinard [26] and Vinard et al. [38] report an underpressured marl-shale aquitard at a depth of 900 m in Wellenberg, Switzerland. Vinard et al. [39] suggest the underpressures could be related to stress relief due to deglaciation, extensive erosion, or tectonic-thrusting scenarios that would result in pore dilation. They also suggest that the presence of a gas phase might contribute to fluid underpressures in the aquitard. Through hydromechanical coupling, changes in in situ stresses affect pore fluid pressures. Johnston et al. [25] and Neuzil and Provost [40] examined the incremental stress induced by crustal flexure which occurs when the lithosphere bends beneath glacial loads. Depending on the state of the groundwater system when deglaciation occurs, dilation of the rock matrix following deglaciation may be sufficient to contribute to the generation of pore fluid underpressures.
Neuzil and Pollock [41] and Jiao and Zheng [27] demonstrate that underpressures caused by the erosion of sedimentary basins can persist for geologically significant periods of time. According to Coniglio and Williams-Jones [42], a minimum of 1500 m of Paleozoic strata is estimated to have been present in the Manitoulin Island area (north-eastern portion of the Michigan Basin) based on data from Cercone [43] (central and northern Michigan Basin). Wang et al. [44] provides a similar estimate of 1500 m of erosion, but for the south-central region of the Michigan Basin. Based on these two studies of burial history in the Michigan Basin, a reasonable maximum of 1000 m erosion of Paleozoic sediment has occurred at the Bruce Site over an approximate 20 million year period [4].
1.1. Prior Analyses of Bruce Site Underpressures
A three-dimensional regional scale model was developed for the purpose of the Bruce Site DGR project encompassing an area of approximately 18500 km2 to investigate, in part, impacts of glaciation on the groundwater system [37, 45]. The analyses honoured the horizontal hydraulic conductivities estimated in the DGR field program while storage coefficients were based on laboratory derived mechanical properties of the rocks at the site [4]. Sykes et al. [37] did not investigate scenarios with alternate parameters. The vertical hydraulic conductivities are based on field estimated horizontal hydraulic conductivities and an assumed horizontal : vertical anisotropy ratio of 10 : 1 for most formations. Many formations have horizontal permeabilities on the order of 10−21 m2 (~10−14 m/s). Laboratory testing of DGR site core from the proposed repository horizon by Selvadurai and Jenner [46] and the testing of samples from the same argillaceous formation closer to surface (Lindsay Formation) yield permeability values between 1.0 × 10−22 (~10−15 m/s) and 1.0 × 10−19 m2 (~10−12 m/s).
A three-dimensional model was deemed necessary to allow the Niagaran Group (comprising the Guelph, Goat Island, Gasport, and Lions Head formations; refer to Figure 2) to subcrop within the model domain and to provide a hydraulic connection for subglacial fluid pressures to propagate to depth. The dimensionality and geometry of geologic units were imported from a three-dimensional constrained geologic reconstruction by ITASCA CANADA and AECOM [47], and the surface hydraulic and mechanical boundary conditions were imported from a detailed glaciation reconstruction by Peltier [48]. As such, no allowance was made for adjusting parameters through a calibration process. Ten different paleohydrogeologic numerical models, which included permafrost and density-dependent groundwater and brine transport were examined. The measured underpressures in the Ordovician formations could not be described by the analyses of Sykes et al. [37]; however, underpressures did form in the Silurian and in the Middle Ordovician in some simulations.
Following the work of Sykes et al. [37], the two-dimensional model studies of Nasir et al. [29] and Khader and Novakowski [30] illustrated that formation underpressures could be partially attributed to the removal of glacial ice loads. The authors suggested that discrepancies between the model predicted results and observed underpressures could be the result of factors not evaluated, such as the presence of a gas phase, osmotic pressure, and/or erosion. Neuzil and Provost [40] and Neuzil [20] stated that glacial forcing was a possible cause of the underpressures based on a one-dimensional analysis at the Bruce Site using a simplified lithology; however, other mechanisms could not be ruled out.
This paper describes analyses that include a detailed characterization of the complete lithology (16 layers) between the Guelph Formation and Cambrian units. The Cambrian units, dominated by sandstones, subcrop to the west and north of Lake Michigan and pinch out against the Algonquin Arch to the east of the Bruce Site (the Cambrian is absent in reference well T004854, located approximately 20 km southeast of the Bruce Site) [49, 50]. The Guelph Formation forms an extensive regional confined aquifer which outcrops along the Bruce Peninsula (approximately 60 km northeast of the Bruce Site) and in Guelph, Ontario (approximately 140 km southeast of the Bruce Site) [50]. Uncertainty in glaciation scenarios is accommodated using two reconstructions (nn9921 and nn9930) of the most recent 120 ka glaciation over the North American continent by Peltier [48]. Uncertainty in model parameterization is accommodated by including horizontal : vertical hydraulic conductivity ratios and vertical compressibility for each layer. Because vertical compressibilities are used as calibration parameters, specific storage and loading efficiency are computed quantities and are not treated independently. Uncertainty in the bottom boundary condition is included as a model parameter, whereas the top boundary condition is considered by assigning varying degrees of hydraulic connectivity between the heads in the Guelph Formation and the subglacial fluid pressures at the Bruce Site. A total of 33 model parameters are estimated for each of the 20 paleohydrogeologic simulations. Furthermore, transport simulations of a helium tracer over 260 Ma, including exhumation of 1000 m of Paleozoic sediments over a 20 Ma period, using the calibrated model parameters from the paleohydrogeologic simulations, are used to identify calibrations which can preserve the distinct helium profile and underpressures over Paleozoic time scales.
2. Bruce Site Characteristics
The sedimentary sequence underlying the Bruce Site (Figure 2) can be partitioned into three regimes based on the rock and groundwater properties. (i) A shallow groundwater zone (0–170 mbgs) is typified by high horizontal hydraulic conductivities (10−5-10−6 m/s), fresh to brackish groundwater and hydraulic gradients governed by topography [4]. The sedimentary rock in this shallow groundwater regime is comprised of the dolostone formations of Devonian and upper Silurian age. (ii) An intermediate groundwater regime (170–470 mbgs) is described as an aquitard that behaves as a transition zone between the shallow and deep groundwater regimes. This intermediate regime is comprised of the upper Silurian Salina Group, the Niagaran Group (including the Guelph, Goat Island, Gasport, and Lions Head formations) and the lower Silurian carbonates and shales. With the exception of the Salina A1 member carbonate and the Guelph Formation, the horizontal hydraulic conductivities of these bedrock formations are typically less than 10−11 m/s. The groundwater and pore water chemistry is typified by brines with the highest total dissolved solids concentrations (370 g/L) occurring in the Guelph Formation. (iii) A deep groundwater regime (470–840 mbgs) includes the Ordovician shales and carbonate formations and the thin basal confined Cambrian sandstone aquifer. The horizontal hydraulic conductivities within the Ordovician sediments range between 10−11 and 10−15 m/s. Pore water within the Ordovician sediments is brines with total dissolved solids concentrations ranging between 220 and 300 g/L.
Subsurface characterization at the Bruce Site was accomplished through field and laboratory investigations at three triangularly located vertical DGR-series boreholes (150 mm diameter) intersecting the intermediate and deep groundwater regimes with a separation distance of 1.2 km. In situ borehole testing yielded estimates of horizontal hydraulic conductivities, and the installation of a Westbay MP55 multilevel casing system in the boreholes yielded in situ pore pressure measurements. An example of the time series head measurements, corrected for fluid density to estimate environmental head, as obtained at borehole DGR-4 is shown in Figure 3. The environmental head plot, as described by Lusczynski [51], is calculated from the measured formation pressures by subtracting the excess heads imposed by pore water with a density greater than freshwater. Pore water densities are determined from site-specific pore water TDS versus density relationships. Given a ground surface elevation of 181.6 mASL, the environmental head profile, shown in Figure 3, indicates that the Ordovician formations are significantly underpressured with respect to ground surface. The freshwater head difference between a hydrostatic condition and the measured underpressures is approximately 300 m; the freshwater head in the Cobourg formation is significantly less than that of Lake Huron (176 masl). The Cambrian unit was observed to be overpressured relative to ground surface by approximately 180 m. The relatively high hydraulic gradients (1-2) vertically convergent towards the Ordovician formations (see Figure 3) are indicative of a transient groundwater system that has not achieved hydrostatic equilibrium. The pore water environmental heads at DGR-4 on November 4, 2013, are used to calibrate the one-dimensional numerical models (see Figure 3).
3. Governing Equations
According to Bear [52], for flow in a saturated porous medium, the equation of mass conservation expressed in terms of hydraulic head iswhere is the hydraulic head; is the vertical hydraulic conductivity; is the fluid source/sink term; is the coefficient of specific storage; and is the time. If the solids of the porous media are considered incompressible, the specific storage and loading efficiency can be expressed aswhere is the fluid density; is the acceleration due to gravity; is the porosity of the porous media; is the fluid compressibility; and is the coefficient of vertical compressibility for the porous media, which can be determined from rock mechanics properties such as Young’s modulus , Poisson’s ratio , and the bulk modulus as [24, 53]Therefore, the one-dimensional vertical groundwater flow equation including hydromechanical coupling iswhere is the vertical stress. The hydromechanical term can be regarded as an additional source/sink term depending on the rate of vertical stress change , loading efficiency, and specific storage. By assuming the loading and unloading are laterally homogeneous and strains are limited to the vertical direction, the specific storage and loading efficiency, accounting for the hydromechanical coupling due to glaciation or exhumation, are dependent on the Biot coefficient, the Poisson’s ratio, and the compressibilities of the porous media, solid grains, and the pore fluid [24, 54]. Mechanical loading due to glacial ice-sheet advance and retreat can influence an underlying groundwater system through the coupling of pore water pressure and a deformable rock matrix. Anomalous hydraulic pressures in low permeability formations can occur, for example, as a remnant of mechanical and/or hydrologic perturbations due to glacial episodes [55]. The loading efficiency ranges from zero to unity, representing the fraction of the applied stress change, which results in a pore fluid pressure change. Underpressures can occur when the pore water outflow during glacial loading is greater than inflow during glacial unloading. This equation demonstrates that underpressures can occur when the rate of glacial unloading is greater than the rate of glacial loading.
For helium transport, the generalized solute transport equation for a saturated porous media is [52]where is the hydrodynamic dispersion tensor; is the concentration; is the Darcy flux; and is the porosity. is further defined by Bear [52] aswhere is the longitudinal dispersivity, is the tortuosity, and is the molecular diffusion coefficient.
4. Modelling Approach
The formation of underpressures and helium transport in the Ordovician sediments is investigated using the computational model HydroGeoSphere [56] and a two-dimensional transient finite element groundwater model with hydromechanical coupling. The latter finite element model was validated against HydroGeoSphere and analytical solutions and can be used for the 1D vertical analyses in this paper by ensuring that boundary conditions and model properties are horizontally constant at any point in time. This latter model was developed specifically for this study to execute ten paleohydrogeologic simulations in less than 5 seconds, as compared to approximately 35 minutes with HGS. As over 627,000 model evaluations were performed, the 2D finite element model was necessary to make this study computationally tractable. Inverse modelling used the BOBYQA algorithm [57]. Implicit time weighting was used for the transient simulations and the LAPACK matrix library was used to solve the resulting system of matrix equations [58]. Fully saturated freshwater analyses are performed; the use of a freshwater analysis is justified as the excess head because a fluid density greater than freshwater can be subtracted by computing the environmental head [51].
Referring to Figure 2, the modelling domain extends vertically from the Guelph Formation to the basal Cambrian units and is comprised of 1862 triangular elements, each with a height of approximately 0.5 m. This domain is exclusively comprised of low hydraulic conductivity formations where hydraulic gradients are predominantly vertical, thereby allowing a one-dimensional modelling domain to be used. The top and bottom boundary conditions for the model were prescribed as temporally varying Dirichlet boundary conditions to replicate measured formation pressures in the bounding aquifers and to accommodate the different conceptual models of glaciation. Table 1, adapted from Sykes et al. [37] and based on the borehole investigations and in situ hydraulic testing, lists thickness and hydrogeologic parameters for each formation in the domain.
For helium transport, a free-water diffusion coefficient of 5.8 × 10−9 m/s2 is used [59]. Tortuosity values are computed from Table of Al et al. [10] for an 18O profile and adjusted to account for formation porosities in Table 1. Longitudinal dispersivity values of 0.5 m are used. The calibrated model parameter values for each of the 20 scenarios are subsequently used in 260 Ma exhumation simulations which include helium tracer transport.
Approximately nine episodes of continental scale glaciation have occurred during the past 1 Ma over the Canadian landmass [60]. Although a future glaciation scenario is of interest, a Bayesian approach is applied to examine a range of models for the most recent North American glaciation event, constrained by various long-term observations in sea level, ice core oxygen isotope ratios, maximal extent of glaciation, and continental isostatic rebound, among others. The physical model used for these simulations is the University of Toronto Glacial Systems Model (GSM) [48]. Peltier [48] describes eight (8) models that provide acceptable fits to the totality of observational constraints. Of these eight models, nn9921 and nn9930 are two of the best models based on aggregate misfit, and both include high-resolution permafrost development. Of the two models, nn9921 and nn9930, model nn9930 had less permafrost than nn9921, yet both are selected for the paleohydrogeologic simulations presented in this paper. Permafrost is excluded from all simulations as the top of the model domain is located at a depth of 370 mbgs, much deeper than the maximum permafrost depth in the GSM models. A plot of ice thickness for the nn9921 and nn9930 GSM model outputs at the Bruce Site is shown in Figure 5. The temporal ice thicknesses are used to assign hydraulic boundary conditions for the top and bottom of the modelling domain. Hydraulic boundary conditions are based on a subglacial freshwater head equivalent to a specified percentage of ice-sheet thickness. The ice-sheet thickness varies in time; however, only the subglacial freshwater head representing the hydraulic boundary condition is scaled, while the mechanical load imposed by the ice-sheet onto the geosphere is not scaled. Glaciation simulation nn9921, shown in Figure 5, provides multiple ice-sheet advance/retreat cycles at the DGR site over the 120,000-year glaciation time period. At approximately 22,000 years before present, the ice-sheet reached a simulated maximum thickness at the Bruce Site of over 3 km.
5. Model Development
To characterize parameter and conceptual uncertainty in the site description, a total of 20 numerical models are developed. Uncertainty in ice-sheet evolution was accommodated using two GSM simulations: nn9921 and nn9930. Uncertainty in model initial conditions prior to glaciation was accommodated by using both a computed steady-state hydrostatic head and by using the present day underpressured environmental heads at the site. Uncertainty in the top boundary condition, representing the subcropping Guelph Formation, was accommodated by varying the percentage of the ice-sheet thickness that contributes to excess head through a range of 0% to 100% in steps of 25%. Although the subglacial head is varied as a function of the ice thickness, the ice-sheet mechanical load imposed on the geosphere is not varied. Having two glaciation scenarios, two initial conditions, and 5 differing percentage values for the excess head generated by glacial loading at the top boundary, a total of 20 models result. Uncertainty in the bottom boundary, representing the Cambrian aquifer, was accommodated by using the percentage of the ice-sheet thickness that contributes to excess head as a calibration parameter. The goal of the modelling was to estimate horizontal : vertical hydraulic conductivity anisotropy ratios (using the horizontal hydraulic conductivities determined from in situ borehole testing) and vertical compressibilities through an inverse modelling approach.
To establish initial freshwater head conditions for the models, heads for the Guelph and Cambrian aquifers are specified based on measurements in the DGR-4 borehole (see Figure 3); the Guelph Formation is set to a head of 194.3 mASL, nearly at present day hydrostatic levels, and the Cambrian unit is set to a head of 366.7 mASL. Hydraulic heads for intervening layers are determined by either steady-state simulation or by linearly interpolating present day environmental heads with depth at the site. The top and bottom prescribed head boundary conditions for each 120 ka glacial cycle are set to the heads in the Guelph and Cambrian aquifers plus a scaled ice-sheet thickness expressed as equivalent freshwater head. The scaling factor for ice-sheet thickness applied at the top boundary condition is fixed for each simulation, while the bottom boundary condition is estimated in the calibration process and is associated with excess hydraulic head at the base of the ice-sheet. This scaling factor for hydraulic head can vary depending on conditions at the interface between the geosphere and an ice-sheet and on the hydraulic connectivity between the bounding aquifers at the Bruce Site and the ground surface. The underpressures in the Ordovician may be the result of multiple glacial cycles of 120 ka each. The paleohydrogeologic simulations in this study apply ten consecutive cycles of either the nn9921 or nn9930 glacial models. The initial condition for each paleohydrogeologic cycle uses the final freshwater heads from the previous paleohydrogeologic cycle.
For each of the 20 models, the horizontal : vertical hydraulic conductivity anisotropy ratios and the vertical compressibility are estimated for each of the 16 layers in the model domain between the Guelph and Cambrian aquifers, for a total of 32 lithologic parameters. The excess head as a percentage of the applied ice-sheet load for the Cambrian aquifer represents an additional calibration parameter. The vertical hydraulic conductivity is bounded based on the measured horizontal hydraulic conductivity and a horizontal : vertical (H : V) anisotropy ratio that varies from 100,000 : 1 to 1 : 1. The specific storage and loading efficiency are calculated from the vertical compressibility . The vertical compressibility is bounded from 10−8 to 10−14 to ensure a computed loading efficiency in a range between 0 and 1. The bulk fluid modulus follows Sykes et al. [37] and is set to 3.3 GPa. The starting point for all calibrations was at the midpoint of the parameter range. The BOBYQA (Bounded Optimization BY Quadratic Approximation) algorithm was chosen to estimate 33 parameters for each of the 20 model simulations. BOBYQA allows for constrained optimization without the need to compute derivatives leading to more stable operation and can accommodate large numbers of parameters [57]. The calibration algorithm determined the set of parameters which minimized the sum of the squared differences between the simulated and observed heads measured in the DGR-4 borehole on November 3, 2013 (see Figure 3). Due to the significant variability in measured heads in both the Coboconk and Gull River units, head measurements within these units were replaced with soft data points, linearly interpolated between the head values in the Cambrian aquifer and the base of the Kirkfield formation, for use in calibration.
In addition, a 260 Ma exhumation scenario, in which 1000 m of erosion occurs over a 20 Ma period, with a unit tracer representing helium isotope ratios was simulated for each of the 20 paleohydrogeologic calibrations with tortuosity values computed from Al et al. [10]. The purpose of these simulations is to determine whether the calibrated parameters from a glaciation simulation (i) can lead to the formation of residual underpressures at the end of the 260 Ma simulation and (ii) whether the helium tracer profile will remain stagnant over the 260 Ma.
6. Simulation Results and Discussion
The evolution of the Ordovician sediments over geologic time has included burial, exhumation, and repeated episodes of glaciation by continental scale ice-sheets [42, 44, 60]. Following exhumation, whereby an estimated maximum of 1000 m of Paleozoic sediment was eroded at the Bruce Site over an approximate 20-million-year period prior to the mid-Jurassic [4], Pleistocene glaciation imposed significant hydraulic and mechanical loads on the geosphere leading to the generation of both over- and underpressures, depending on the rock mechanical properties and the rate of glacial loading and unloading [20, 29, 37, 40].
Computed heads at the end of ten paleohydrogeologic simulations for glaciation scenario nn9921 and varying percentages of excess head for the top boundary condition are shown in Figure 6. Figure 6(a) uses present day environmental heads as the initial condition while Figure 6(b) uses a hydrostatic head condition computed assuming a steady-state groundwater system. Figure 7 shows final heads using the nn9930 glaciation simulation. The calibrated values of horizontal : vertical hydraulic conductivity anisotropy ratios for glaciation scenarios nn9921 and nn9930 are shown in Figures S1 and S2, respectively, while the calibrated values of vertical compressibility for glaciation scenarios nn9921 and nn9930 are shown in Figures S3 and S4, respectively. The impact of both the anisotropy ratios and the vertical compressibility can be more readily assessed by considering the computed values of vertical hydraulic conductivity (Figures 8 and 9), specific storage (Figures S5 and S6), and loading efficiency (Figures S7 and S8). The calibrated values for the Cambrian aquifer boundary condition ranged from 13% to 25% (of ice-sheet thickness expressed as equivalent freshwater head) for six of the nn9921 based paleohydrogeologic simulations and 0% for the remaining four nn9921 simulations. For the nn9930 based paleohydrogeologic simulations, all calibrated values for the Cambrian aquifer boundary condition were 0% of ice-sheet thickness expressed as equivalent freshwater head.
Based on the calibration results, it can be concluded that the observed underpressures can be described by mechanical loading during glaciation and unloading and consequent rock dilation during deglaciation; however, the calibrated parameters depend on the choice of glaciation boundary conditions, initial head conditions, and the amount of excess head from glaciation that is applied at the top and bottom boundaries of the models. Once underpressured, the pattern of underpressures after each successive glacial cycle was repeated and may indicate that an underpressured condition preceding glaciation, caused by exhumation, gas generation, or some other mechanism, can be maintained through multiple glacial cycles.
In comparing the horizontal : vertical hydraulic conductivity anisotropy ratios for the nn9921 (Figure S1) and nn9930 (Figure S2) glaciation simulations, the simulations with initial conditions representing present day underpressures require lower vertical hydraulic conductivities in the upper Ordovician than for an initial hydrostatic condition. Lower vertical hydraulic conductivities, by up to 5 orders-of-magnitude (see Table 2), are needed to maintain the system in an underpressured state regardless of the glaciation simulation used, whereas higher vertical hydraulic conductivities are needed to allow pore water in the Georgian Bay and Blue Mountain formations to migrate out of the groundwater system from the initial hydrostatic condition towards an underpressured condition. The middle Ordovician and upper Ordovician are distinctly different in their response to glaciation, likely attributed to pore fluids exiting or entering the groundwater system either from the top or from the bottom boundaries. This can also be seen in the calibrated vertical compressibilities which are generally more tightly grouped in the upper Ordovician as compared to the middle Ordovician.
The impact on computed heads of the initial condition can be seen in Figures 6 and 7 for the nn9921 and nn9930 glaciation simulations, respectively. For greater numbers of glacial advances and retreats (nn9921, Figure 6), the hydrostatic initial condition simulation is able to reasonably approximate the present day head profile at the Bruce Site for the 0%–75% top boundary simulations. Only the 100% simulation is significantly different. In the case of the nn9930 glaciation scenario, only the 0%–25% top boundary scenarios are able to reasonably approximate the present day head profile at the Bruce Site, primarily due to perfect sink created at the top boundary by ignoring excess head from glaciation. This demonstrates that the calibrated system properties are more sensitive to glaciation reconstructions used to assign the paleohydrogeologic boundary conditions if the groundwater system is initially in a hydrostatic state. The robustness and resilience of the groundwater system to external perturbations is greater for the state whereby underpressured conditions predate the onset of glaciation.
The key requirement to replicate the underpressures in the analyses is that the Guelph and Cambrian head boundary conditions can only represent a fraction of the ice thickness to allow drainage to occur. Under such conditions, hydraulic gradients leading to pore fluid exiting the groundwater system can be generated during glacial loading through hydromechanical coupling, which contributes to the development of underpressures at the end of paleohydrogeologic simulation cycle. For the bottom boundary condition, the applied freshwater head for four of the nn9921 based simulations and all of the nn9930 based simulations is 0% of the ice-sheet load or is essentially independent of the ice-sheet loading and less than 25% for the remaining six scenarios. In such conditions, the Cambrian aquifer acts as a drain for in situ generated pore pressures. The Cambrian units have a very high horizontal hydraulic conductivity of 3.0 × 10−6 m/s [5, 61] which can dissipate transient hydraulic perturbations efficiently. Similar findings were reported by Khader and Novakowski [30] using a two-dimensional model. The analyses of Sykes et al. [37] included one simulation that allowed the Cambrian aquifer to drain and an underpressure in the middle Ordovician formed by the end of one glacial cycle; however, additional glacial cycles were not simulated.
Clark et al. [9] present hydrogeochemical and isotopic pore water data for helium (see Figure 4) within the Ordovician shales and carbonates that reveal an essentially immobile groundwater system over a period of hundreds of millions of years. A shift in 3He/4He ratios occurs at the base of the halite-mineralized section of the Ordovician shales and Cobourg limestone, further indication of a stagnant and immobile groundwater system [9]. The presence of halite or a second fluid phase (horizontal gas lenses or lamina of oil) identified at the site [4] can reduce effective vertical hydraulic conductivities through porosity reduction [9] or multiphase flow [62], respectively. Al et al. [10] describe the geochemical evolution of the Bruce nuclear site and potential mechanisms for halite precipitation, including hyperfiltration and temperature decreases during consolidation of the Ordovician shales following maximum burial. With the aid of numerical modelling of Cl and 18O tracers, Al et al. [10] suggest that effective diffusion coefficients for the Ordovician carbonates are up to an order-of-magnitude lower than laboratory estimated values and would act to further decrease vertical migration. As pore fluid fluxes are vertically oriented within the horizontally bedded sediments, a significant reduction in the permeability of even a relatively thin section of the domain can influence the entire domain due to the consequence of harmonic averaging [63].
The simulated tracer profile over a period of 260 Ma did not preserve the step change in the observed helium isotope profile using the tortuosity values in Table 1. As such the tortuosity of the Cobourg and Sherman-Fall formations and hence its effective diffusion coefficient was lowered by 3 orders-of-magnitude in order to preserve a step change in the helium tracer profile over 260 Ma. This approach was considered reasonable given the presence of an aquiclude zone as identified by Clark et al. [9]. Figures 10–13 summarize the simulation of the tracer over 260 Ma using the calibrated vertical hydraulic conductivities from the paleohydrogeologic simulations.
In comparing the figures, the simulations with initial underpressures (Figures 10 and 12) demonstrate less variability in tracer migration as compared to the initial hydrostatic pressure simulations (Figures 11 and 13). This variability is due to the greater advection of the tracer resulting from increased vertical hydraulic conductivities in the simulations which are initially at hydrostatic pressure conditions. These simulations also generate underpressures at the end of the 260 Ma, although the underpressures do not match the present day environmental head profile. A closer match to the field measured underpressures would be obtained by increasing the exhumation interval from its present value of 20 Ma. In addition, regardless of the reduced tortuosity, the nn9930 simulations with the 25% glacial top boundary condition and hydrostatic initial head conditions exhibit sufficient advective transport to eliminate most of the step change in the simulated tracer profile (see Figure 13). In order to preserve a step change in the helium profile with depth, both very low vertical hydraulic conductivities and very low effective diffusion coefficients are needed. The former suggests that underpressures predate glaciation, possibly generated through erosion of up to 1000 m of Paleozoic sediments. The simulations in this study demonstrate that preglacial underpressures, caused by an exhumation event in the Paleozoic, can be preserved through multiple glacial cycles to the present. The simulations in which the initial heads are underpressured result in better preserving the present day helium profile.
7. Conclusions
The primary goal of the study was to explore, through a series of numerical simulations constrained by detailed field observation, factors affecting the nature of long-term phenomena responsible for the generation and preservation of formation underpressures. Three families of numerical experiments for underpressure formation were examined by means of one-dimensional hydromechanically coupled models through the hydrostratigraphic column: (i) uncertainty in glaciation scenarios; (ii) uncertainty in initial heads; and (iii) uncertainty in the degree of hydraulic connectivity between the Guelph Formation at the Bruce Site and the applied glacial loading. The paleohydrogeologic simulations assumed fully saturated conditions. This study aimed to address these uncertainties through qualitatively comparing the inverse model results.
For the bottom boundary condition, the applied freshwater head for most simulations was 0% of the ice-sheet load. As such, the role of the Cambrian aquifer is essentially independent of glacial loading as the formation must act as a drain for in situ generated pore pressures, whether the groundwater system is initially underpressured or hydrostatic. The Cambrian units have a very high horizontal hydraulic conductivity of 3.0 × 10−6 m/s [5, 61] which can dissipate transient hydraulic perturbations efficiently; however, it pinches out within a distance less than 20 km east of the site and it subcrops west of Lake Michigan [49].
Of the two glaciation simulations (nn9921 and nn9930), nn9921 was better able to reproduce the present day pore pressure profile at the DGR site, across a range of glacial boundary conditions imposed at the top boundary and initial conditions. Underpressured initial heads for the paleohydrogeologic simulations lead to lower vertical hydraulic conductivities. All of the simulations resulted in calibrated vertical hydraulic conductivities that were essentially at or below the horizontal hydraulic conductivities from field and laboratory testing at the Bruce Site and applied in the paleohydrogeologic modelling of Sykes et al. [37]. All of the simulations with underpressured initial heads resulted in good fits to the observed head profile, whereas most of simulations with hydrostatic initial heads yielded poorer fits to observed heads.
Evidence exists for very low hydraulic conductivities in specific horizons of the hydrostratigraphic column. For example, the immobility of some environmental tracers such as helium, as documented by Clark et al. [9], demonstrates the existence of very low vertical hydraulic conductivity (≪10−14 m/s) zones that also exhibit very low effective diffusion coefficients. The use of site-specific analogues (self-analogues such as helium) is important for characterizing the future behaviour expected to cause the transport of radionuclides from the site [21]. Due to harmonic averaging, relatively small, but contiguous, layers of reduced hydraulic conductivity can significantly lower effective vertical hydraulic conductivities and dominate the behaviour of vertical pore fluid migration. It should be noted that low vertical hydraulic conductivities (<10−12 per Neuzil [20]) are necessary in order to generate and preserve underpressures [3, 36]. Exhumation analyses were conducted whereby a helium tracer was simulated for each of the 20 calibrated paleohydrogeologic models over a period of 260 Ma. Simulations with initial underpressures (Figures 10 and 12) demonstrate less variability in tracer migration and less advection, as compared to the initial hydrostatic pressure simulations (Figures 11 and 13). These simulations also result in underpressured conditions at the end of the 260 Ma period indicating that the lower vertical hydraulic conductivity values resulting from the initially underpressured paleohydrogeologic scenarios are compatible with generating underpressures through exhumation. The degree of underpressure can be reduced to better match measured values by increasing the exhumation interval from its present value of 20 Ma.
Although glaciation has been demonstrated to generate underpressures in the literature, this study demonstrates that other mechanisms, such as exhumation, can generate underpressures and that these underpressures, predating glaciation, can be preserved through repeated cycles of glaciation. In order for this to occur, vertical hydraulic conductivities need to be lower than if glaciation were considered as the sole mechanism. The lower vertical hydraulic conductivities are also needed to better preserve a step change in the simulated helium tracer profiles across all scenarios over a time scale of 260 Ma. The very low vertical hydraulic conductivities necessary for generating and preserving underpressures through exhumation over a period of 260 Ma, preserving a helium tracer profile over that same time, and preserving underpressures through ten glacial cycles, together suggest that underpressures predate glaciation. The robustness and resilience of the groundwater system to external perturbations are greater for the state where underpressured conditions predate the onset of glaciation. Although these analyses were focused on the Bruce Site, the methodology and findings are applicable to any site that has undergone multiple mechanisms in the formation of in situ underpressures, such as exhumation and repeated cycles of glaciation.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this article.
Acknowledgments
This study was funded by the Nuclear Waste Management Organization (NWMO) of Canada. Data for the models developed herein are provided in the included tables.
Supplementary Materials
Figure S1: calibrated horizontal : vertical hydraulic conductivity anisotropy ratios using glaciation scenario nn9921. Coloured lines represent different percentages of ice-sheet thickness applied as excess head to the top boundary condition. Subfigure (a) uses present day environmental heads as the initial condition while subfigure (b) uses a hydrostatic head condition computed assuming steady-state. Figure S2: calibrated horizontal : vertical hydraulic conductivity anisotropy ratios using glaciation scenario nn9930. Coloured lines represent different percentages of ice-sheet thickness applied as excess head to the top boundary condition. Subfigure (a) uses present day environmental heads as the initial condition while subfigure (b) uses a hydrostatic head condition computed assuming steady-state. Figure S3: calibrated vertical compressibility using glaciation scenario nn9921. Coloured lines represent different percentages of ice-sheet thickness applied as excess head to the top boundary condition. Subfigure (a) uses present day environmental heads as the initial condition while subfigure (b) uses a hydrostatic head condition computed assuming steady-state. Figure S4: calibrated vertical compressibility using glaciation scenario nn9930. Coloured lines represent different percentages of ice-sheet thickness applied as excess head to the top boundary condition. Subfigure (a) uses present day environmental heads as the initial condition while subfigure (b) uses a hydrostatic head condition computed assuming steady-state. Figure S5: calibrated specific storage using glaciation scenario nn9921. Coloured lines represent different percentages of ice-sheet thickness applied as excess head to the top boundary condition. Subfigure (a) uses present day environmental heads as the initial condition while subfigure (b) uses a hydrostatic head condition computed assuming steady-state. Figure S6: calibrated specific storage using glaciation scenario nn9930. Coloured lines represent different percentages of ice-sheet thickness applied as excess head to the top boundary condition. Subfigure (a) uses present day environmental heads as the initial condition while subfigure (b) uses a hydrostatic head condition computed assuming steady-state. Figure S7: calibrated loading efficiency using glaciation scenario nn9921. Coloured lines represent different percentages of ice-sheet thickness applied as excess head to the top boundary condition. Subfigure (a) uses present day environmental heads as the initial condition while subfigure (b) uses a hydrostatic head condition computed assuming steady-state. Figure S8: calibrated loading efficiency using glaciation scenario nn9930. Coloured lines represent different percentages of ice-sheet thickness applied as excess head to the top boundary condition. Subfigure (a) uses present day environmental heads as the initial condition while subfigure (b) uses a hydrostatic head condition computed assuming steady-state. (Supplementary Materials)