Representing Glaciations and Subglacial Processes in Hydrogeological Models: A Numerical Investigation

Sterckx, Arnaud; Lemieux, Jean-Michel; Vaikmäe, Rein

doi:https://doi.org/10.1155/2017/4598902

Geofluids

On this page

Abstract Introduction Methods Results Discussion Conclusion Disclosure Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2017 | Article ID 4598902 | https://doi.org/10.1155/2017/4598902

Representing Glaciations and Subglacial Processes in Hydrogeological Models: A Numerical Investigation

Arnaud Sterckx,^1,2Jean-Michel Lemieux,^1,2and Rein Vaikmäe³

Academic Editor: Marco Petitta

Received16 Dec 2016

Accepted28 Feb 2017

Published04 Apr 2017

Abstract

The specific impact of glacial processes on groundwater flow and solute transport under ice-sheets was determined by means of numerical simulations. Groundwater flow and the transport of δ¹⁸O, TDS, and groundwater age were simulated in a generic sedimentary basin during a single glacial event followed by a postglacial period. Results show that simulating subglacial recharge with a fixed flux boundary condition is relevant only for small fluxes, which could be the case under partially wet-based ice-sheets. Glacial loading decreases overpressures, which appear only in thick and low hydraulic diffusivity layers. If subglacial recharge is low, glacial loading can lead to underpressures after the retreat of the ice-sheet. Isostasy reduces considerably the infiltration of meltwater and the groundwater flow rates. Below permafrost, groundwater flow is reduced under the ice-sheet but is enhanced beyond the ice-sheet front. Accounting for salinity-dependent density reduces the infiltration of meltwater at depth. This study shows that each glacial process is potentially relevant in models of subglacial groundwater flow and solute transport. It provides a good basis for building and interpreting such models in the future.

1. Introduction

Glaciations are known to have a large and long-lasting impact on groundwater flow. Several geochemical and isotopic studies have shown that water of glacial origin is still present in basins in North America (e.g., [1, 2]) and northern Europe [3] that were formerly covered by ice-sheets (see [4] for a review of these studies). The presence of glacial water is commonly identified by low salinity, depletion in heavy isotopes such as ¹⁸O and ²H, high excess air, cold recharge temperature inferred from noble gases, and old groundwater age.

Subglacial recharge of meltwater is the prevailing hypothesis of recharge of water of glacial origin during the Pleistocene (e.g., [3, 5]). Ice-sheets can be partially wet-based [6] due to basal friction and can release large volumes of subglacial meltwater. If the permeability of the ground is sufficient, it was put forward that meltwater could infiltrate into the subsurface under the pressure exerted by the ice-sheet [7, 8]. Grundl et al. [9] showed that recharge of supraglacial meltwater can also occur below ice-sheets due to the presence of crevasses. In this case, the pressure at the base of the ice-sheet would not be induced by the weight of the ice but by the weight of the column of water within the ice. The presence of fresh water at depth [2] or in abnormal position with respect to present-day hydraulic gradient [10] suggests that the trend and the intensity of regional flow were considerably modified. Modelling studies confirmed that the injection of large volumes of meltwater has the potential to reorganize entirely the pattern of regional flow [7] and even to reverse the direction of the flow if the ice-sheet runs opposite to the direction of interglacial topography-driven flow [11].

Understanding the impact of glaciations on groundwater flow systems is important for several reasons. First, waters of glacial origin are usually of good quality and are, in some locations, an important source of nonrenewable freshwater [3, 12]. The sound management of these resources requires a good scientific understanding of their origin, their extent, and the driving forces acting on them. Secondly, groundwater flow can be an important transport vector for nuclear waste radionuclides. In cases where the containment mechanism of spent nuclear fuel is compromised, groundwater flow could remobilize the radionuclides into the ecosystems and increase the health risks to fauna and humans who depend on it [13]. Given the long time scale required for spent nuclear fuels to decay to safe levels, the suitability of deep geologic repositories (DGR) must be demonstrated over timescales of 100 ky. According to geological history, it is likely that at least one glaciation event will occur during this period [14]. Therefore, DGR stability and the risks of contamination from groundwater flow during a glaciation event depend on the dynamics of groundwater flow below ice-sheets [15]. Another motivation for studying glaciations in hydrogeology is the migration or the enhancement of hydrocarbon resources in deep basins affected by ice-sheets [16, 17].

In North America and northern Europe, ice-sheets have retreated at the end of the Pleistocene, more than 10 ky BP. In these regions, groundwater flow systems containing glacial water are not at equilibrium with present-day conditions and their study requires a transient approach, at the scale of climate changes [18]. Numerical models of groundwater flow are powerful tools that can be used to quantify Pleistocene recharge origin, rates, and patterns, along with groundwater dynamics. Recent advances in modelling have allowed taking into account the numerous processes involved during a glacial period. The state-of-the-art models now include sophisticated transient 3D groundwater flow, subglacial recharge, permafrost evolution, glacial isostasy, ice loading, and sea-level change [19–27].

However, the paucity of field data as well as uncertainties in paleoclimatic forcing, palaeoboundary conditions, and petrophysical properties pose a major challenge for calibrating and validating these models [15]. Moreover, the amount of processes considered makes it difficult to understand their respective impacts on groundwater flow dynamics. It is also unclear whether all these processes are important to be considered for adequately capturing the impact of glaciations on groundwater flow dynamics.

The objective of this study is to assess the specific impact of glacial processes on groundwater flow dynamics and to find out which are relevant to include into a numerical model for efficiently capturing the impact of glaciations on groundwater flow at the basin scale. To do so, glacial processes that are characteristic of subglacial environments and could have an impact on groundwater flow are first identified and the strategy for implementing those processes into a numerical model is discussed. Those glacial processes are subglacial recharge, poroelastic deformation of rocks under the weight of the ice-sheet, isostasy, surface drainage evolution, permafrost evolution, and density-dependent flow involving glacial meltwater and deep brines. The individual impacts of these processes on groundwater flow dynamics and solute transport are then assessed using a simple conceptual model representing a sedimentary basin. The finding of this study should help determining which of these processes need to be simulated for adequately representing glaciations in hydrogeological numerical models.

2. Glacial Processes

Over the last decades, a wide range of numerical models of subglacial hydrogeology have been published, spanning different space and time scales, being two- or three-dimensional and steady-state or transient, solving solute transport, heat transport, or deformation equations in addition to the flow equation, and, last but not least, representing different glacial processes. This section briefly presents those glacial processes and the challenges associated with their numerical representation.

2.1. Subglacial Recharge

Numerical models have demonstrated that subglacial recharge is able to explain the presence of large volumes of fresh glacial water at great depth or in abnormal position with respect to present-day flow conditions [5]. Therefore, the fundamental relevance of this process on subglacial hydrogeology is not discussed in this paper. What is discussed instead is the impact of the boundary condition used for representing this process numerically. Two approaches are commonly used: a prescribed head below the ice-sheet or a specified flux representing the inflow of meltwater into the subsurface.

When a fixed head is chosen, it is most of the time considered that the head is equivalent to the weight of the ice-sheet, after the following relation:where is the head prescribed beneath the ice-sheet [m], 0.9 is the ratio between the density of ice and the density of liquid water, and is the ice-sheet thickness [m] (e.g., [19, 21, 22, 26, 27]). Whether this head represents the weight of the ice-sheet or the water column within the ice, this boundary condition can be seen as an upper limit, because any upper value would cause the ice-sheet to “float.” For this reason, some authors [23, 28] have assigned smaller heads, which means that not all the weight of the ice-sheet was supported by basal meltwater or the water-table was located lower within the ice-sheet. However, the rare field measurements of subglacial pressure under modern ice-sheets support the use of (1) [29, 30].

When a specified flux is used, it can be constrained by a maximal head equivalent to the weight of the ice [20, 24]: fixed flux turns into fixed head when basal meltwater reaches the pressure of the ice-sheet. In the literature, meltwater fluxes usually vary between 0 and 30 mm/yr (e.g., [5, 11, 32]).

2.2. Glacial Loading

One of the glacial processes much subject to speculation is the deformation of the porous medium. Neuzil [33] identified two types of deformation: direct loading (or compaction) under the weight of the ice and flexural loading in subglacial and periglacial regions affected by glacioisostatic flexure of the lithosphere. Direct loading is supposed to generate abnormal pressures in low hydraulic diffusivity porous media, where pressure changes induced by direct loading/unloading dissipate slowly [19]. Nasir et al. [34], Khader and Novakowski [25], and Neuzil and Provost [35] showed that direct loading could explain abnormal underpressures observed in the eastern Michigan Basin, USA. The effect of flexural loading is more obscure, because flexure causes compressive and extensive stresses whose distribution and intensity evolve transiently. Neuzil [33] and Khader and Novakowski [25] argue that flexural loading could be of greater magnitude than direct loading. However, no clear relation between abnormal pressures and flexural loading was found on the field [33]. Given the complexity to represent flexural loading numerically, few hydrogeological models have considered this process [25].

For this reason, only direct loading was considered in the present study. To simulate direct loading, constitutive equations for poroelastic deformation and flow should be coupled, which can be done numerically [36]. However, in a review of the hydromechanical effects of glaciations, Neuzil [33] notes that the poroelastic deformation equation requires mechanical properties as well as initial and boundary conditions that are poorly known. In fact, most of the studies used a simplified, uncoupled approach in which the deformation of the porous medium is only vertical [19–21, 23, 26, 35]. With this method, any increment of ice-sheet thickness produces an increment of hydraulic head. Other than assuming elastic deformation, this approach supposes that the deformation is linear. Linearity requires that the stresses are small and that the variations of compressibility and permeability can be neglected. According to Neuzil [33], linear elasticity is an acceptable assumption if the porous medium has experienced preconsolidation during previous glaciations. Linear elasticity could be unfit for describing the deformation of fractured media.

2.3. Isostasy

Isostasy refers to the gravitational equilibrium between the lithosphere and the underlying mantle. The current isostatic uplift following the Last Glacial Maximum is generally well documented by paleoshorelines. However, evidences of isostatic depletion during and prior to the Last Glacial Maximum are sparse. Past topographic elevations estimated by inverse models can be used, if available [20, 25]. Otherwise, isostatic depletion and rebound can be assumed to be proportional to the ice-sheet thickness [19, 21, 26]. Considering several kilometers thick ice-sheets, the change in surface elevation could reach several hundreds of meters in places. This change being spatially heterogeneous, the regional hydraulic gradients are expected to be considerably modified [20].

2.4. Permafrost

Permafrost is a soil whose temperature remains below 0°C for at least 2 consecutive years. With some exceptions, it means that at least a fraction of the water content is frozen. The presence of frozen water decreases the effective porosity of the medium and its permeability [37, 38]. Numerical simulations having considered permafrost have shown that it may behave as a confining unit and promote subpermafrost overpressure development in the vicinity of the front of the ice-sheet (e.g., [5, 11, 32, 39]). For this reason, permafrost has been considered in several studies (e.g., [19, 20, 23, 24, 27]).

2.5. Surface Drainage

Large bodies of surface water (rivers, seas, and lake) are of importance in any groundwater flow system. In the context of glaciations, they can be particularly relevant for at least two other reasons. First, in permafrost areas, the bed of sufficiently thick water bodies can be preserved from freezing and promote talik development, ensuring a pathway for recharging or discharging groundwater [7, 32]. Secondly, during ice-sheet melting, large proglacial water bodies are usual because large volumes of meltwater are released that fill in the isostatic depression. During the period following the retreat of the ice-sheet, the level of these proglacial water bodies generally decreases in relation to the postglacial isostatic uplift. Since the rate of uplift varies spatially, the history of shorelines displacement varies from one location to another (e.g., [40]). In places, this history is further complicated by the rapid drainage of lakes or by marine transgressions. Examples can be found, for instance, in Björck [41] for the Baltic Sea history or in Clarke et al. [42] for the final stage of the Agassiz Lake in North America. It results that the base level of groundwater flow constantly changes, both spatially and temporally. In North America, several models have demonstrated the role of subglacial recharge and sea-level changes in the emplacement of glacial freshwater offshore on the Atlantic Continental Shelf, USA [12, 21, 26]. Similar simulations have also been applied to Greenland [43]. In the present study, the relevance of base level changes in relation to postglacial uplift is considered.

2.6. Variable Density and Viscosity

Many authors have simulated the variation of density and viscosity, because meltwater can be driven to great depth in sedimentary basins where brines are present [19, 20, 22–24, 27] or be in contact with seawater in coastal areas [12, 21, 26, 43]. It was showed that more saline water restrains the emplacement of fresh meltwater into the deeper part of sedimentary basins and into coastal aquifers. In the present study, effects of density-driven flow are assessed. Variations of viscosity however are not considered, for simplicity.

3. Methods

The selected glacial processes were implemented in a transient groundwater flow and solute transport numerical model, representing a generic sedimentary basin. This model is intended to portray subglacial groundwater flow dynamics at the regional scale, where subglacial recharge is the dominant process. It is thus not representative of periglacial regions, where other processes such as permafrost or lithospheric flexure could be more significant. Reactive transport is not considered, as it was already studied by Bea et al. [44].

To assess the respective impact of the selected glacial processes on groundwater flow dynamics, the following procedure was followed. First, a reference case was simulated, considering only subglacial recharge. Then, the other glacial processes were added individually to the reference simulation. The results of each simulation were compared to the reference simulation in order to determine the relative impact of each individual glacial process. Table 1 summarizes all the simulations that were performed in this study.

3.1. Conceptual Model

The conceptual model used for the simulations is a two-dimensional cross-section featuring a generic sedimentary basin (Figure 1). Its geometry is typical of intracratonic sedimentary basins, so that the results obtained here can be useful to understand groundwater flow dynamics at the regional scale in real natural systems. The basin is 700 km long and has a maximum thickness of 2000 m. It lies on an impervious crystalline basement. On the upper boundary, hydraulic heads equivalent to the surface elevation are assigned. Since the topography is assumed to be flat, groundwater is in hydrostatic conditions at the beginning of the simulations. The purpose of this choice is to capture the long-term effects of glaciations (e.g., the dissipation of abnormal water pressures), without interference of case-specific topography-induced regional flow. The sedimentary basin contains three aquifer units (coarse grained sandstones) separated by two aquitard units (shales). Quaternary deposits with intermediate hydraulic conductivity are represented with a 25 m thick layer immediately below the surface. Hydraulic and transport properties of the 6 hydrogeological units are summarized in Table 2. Longitudinal and transversal dispersivity values are 10 and 1 m, respectively. Vertical hydraulic conductivity is 10 times lower than horizontal hydraulic conductivity.

3.2. Reference Simulation

Reference simulation A considers only subglacial recharge. Subglacial recharge is simulated with a fixed head boundary condition, stating for the weight of the ice-sheet, because it represents an upper boundary case of the flow regime at the base of an ice-sheet. Applying a fixed head boundary condition requires knowing the height of the ice-sheet. Perfectly plastic ice-sheets have a semiparabolic profile [45] that is described by a function such aswhere is the ice-sheet thickness [m] at radius [m], is the maximal ice-sheet thickness [m], and is the radius of the front of the ice-sheet [m]. However, the hydraulic head profile is not likely to mirror such a profile, with a sharp gradient at the front of the ice-sheet. In order to smooth the profile of hydraulic heads near the front of the ice-sheet, the following relation was used instead:where is the hydraulic head assigned beneath the ice-sheet [m]. The glaciation scenario considered for this reference simulation is taken from Bense and Person [19] and displayed on Figure 2. At the glacial maximum, half of the basin is covered by the ice-sheet, which reaches a maximal thickness of 1700 m.

Three tracers are used to track meltwater migration in the subsurface: δ¹⁸O, salinity, and groundwater age, which implies solving solute transport equations. The initial and the boundary conditions used for these three solute transport processes are summarized in Table 3.

3.3. Simulations of Individual Glacial Processes

3.3.1. Simulations B1 and B2: Subglacial Recharge Boundary Condition

These two simulations investigate the effect of representing subglacial recharge by a fixed flux boundary condition constrained by the weight of the ice-sheet, instead of assigning a fixed head boundary condition as in reference simulation A. Simulations B1 and B2 assume fluxes of 2 and 10 mm/y, respectively.

3.3.2. Simulation C: Glacial Loading

As in previous studies, the uncoupled approach is followed for simulating glacial loading. Changes in hydraulic heads are proportional to changes in ice-sheet thickness, after the relationwhere is the hydraulic conductivity [m/s], is the pressure of the fluid [Pa], is the water density [kg/m³], is the gravity acceleration [m/s²], is the elevation [m], is the specific storage [m⁻¹], is the loading efficiency [−], and is the pressure exerted by the ice-sheet [Pa]. The loading efficiency [−] is the ratio of the ice-sheet weight transferred to the fluid under undrained conditions, given by the equation [20]where and are the compressibility of the porous medium and the water, respectively [Pa⁻¹]. Here, it is assumed that , which means that 100% of the weight of the ice-sheet is transmitted to groundwater. It represents an upper boundary case of the effect of glacial loading on subglacial groundwater flow.

3.3.3. Simulations D1 and D2: Isostasy

The impact of isostasy is assessed in simulations D1 and D2. Simulation D2 further considers that the postglacial isostatic depression is filled up by a large water body (a lake or a sea) whose level is equal to present-day sea level. During the growth of the ice-sheet, the elevations of the nodes and the hydraulic heads are lowered by 25% of the ice-sheet thickness, similarly to Bense and Person [19]. When the ice-sheet retreats, the uplift is given by the following geometrical relation:where and are the isostatic depletions at two consecutive time steps [m] and is a unit-less factor depending on the length of the time step, so that the uplift has also a postglacial phase. Using for 1 ky time steps yields a maximal postglacial depletion of 80 m at the end of the glaciation and 3 m at the end of the simulation.

3.3.4. Simulation E: Permafrost

In this simulation, the soil is assumed to be initially frozen over 200 m below the surface and over 300 km from the northern end of the basin. When overridden by the ice-sheet, the frozen soil thaws. The following analytical solution from Lunardini [31] is used to calculate the gradual thawing of permafrost:where and are the top and bottom depths of thawed permafrost [m], is the bulk thermal conductivity of thawed soil [W/m·K], is the increase in temperature [K] at the surface causing the degradation of permafrost, is the elapsed time since the increase of temperature [s], and is the latent heat of melting water [J/g]. Basically, these equations consider that the permafrost layer is at 0°C and starts thawing at the top because the surface temperature increases instantaneously by . Simultaneously, the supply of geothermal heat causes the thawing of permafrost from the bottom. The values assigned to the parameters of the solution are listed in Table 4. Since several studies have suggested that the temperature of subglacial recharge was slightly superior to 0°C [3, 9], an initial increase of temperature of 3°C is used. Elements comprised within the permafrost thickness are deactivated, which means that permafrost acts as an impervious layer. No taliks are simulated. Geochemical processes such as isotopic fractionation [46] and salt rejection [47] are not considered. Again, it is a simplified representation of permafrost. It illustrates an upper boundary case of the permafrost impact on groundwater flow.

3.3.5. Simulation F: Variable Density

In this simulation, water density varies with salinity, which promotes vertical buoyancy flow. The effects of temperature and pressure on density are neglected. Water density is assumed to vary linearly with salinity. Measured values of density and total dissolved solids in the Michigan Basin brines reported in McIntosh et al. [22] suggest a relation likewhere and are the densities of saline and fresh water, respectively, [kg/m³] and TDS is the total dissolved solids [kg/m³]. The maximal salinity being 300 g/L in the lower part of the basin, the maximal density of water is 1200 g/L. Accordingly, the initial hydrostatic hydraulic head profile varies from 0 to 200 m in the deeper part of the basin. The vertical salinity gradient can be seen as a compromise between the gradients observed in the Illinois, Michigan, and Appalachian sedimentary basins in North America [48].

3.4. Numerical Implementation

The numerical model FEFLOW [49] was used to conduct these simulations. The code is based on the finite element method and is able to solve 2D or 3D transient problems involving density-dependent groundwater flow and advective/dispersive solute transport. The plug-in interface of FEFLOW was used to implement the glacial processes, for example, to assign boundary conditions and material properties that vary both in time and in space. A triangular mesh of about 300,000 elements was generated. Minimum and maximum surfaces of elements are 20 and 20000 m², respectively. A fully implicit time scheme was used, with 50 y long time steps.

4. Results

In this section, the results of reference simulation A are first presented, followed by the results of the other simulations, where the differences compared to simulation A are highlighted. Figure 3 summarizes the evolution of hydraulic head, δ¹⁸O, salinity, and groundwater age in all simulations, in the two deepest aquifers. A particular attention was given to simulated groundwater flow rates and meltwater volumes, which are shown on Figures 4 and 5.

4.1. Reference Simulation

Figure 6 shows the evolution of hydraulic head, δ¹⁸O, salinity, and groundwater age during simulation A. During the growth of the ice-sheet, hydraulic heads rapidly increase in all layers but the lower aquitard, where the hydraulic diffusivity is the lowest. In this layer, underpressures develop. In overlying aquifers, hydraulic heads mimic the surface conditions, while, in the lowest confined aquifer, hydraulic heads increase beyond the front of the ice-sheet. Due to the logarithmic profile of the ice-sheet, the hydraulic gradient is higher near the front of the ice-sheet, where groundwater flow rates as high as 800 mm/y are observed in the intermediate aquifer (layer 4). Transport of meltwater in the three aquifers is clearly displayed by the contrast in δ¹⁸O, salinity, and groundwater age with respect to pristine groundwater. By the end of the glacial maximum (15 ky BP), the meltwater contents in the intermediate and lower aquifers are about 10 and 12.5%, respectively. In the same aquifers, meltwater has travelled about 50 and 75 km.

During the ice-sheet retreat, hydraulic heads recover to their original static state, except in the lower aquitard, where underpressures turn into overpressures. Maximum overpressures of 800 m are noticed after the final retreat of the ice-sheet. They decrease to 70 m by the end of the simulation. The dissipation of these overpressures leads some meltwater to flow out of the basin but, in the absence of any other flow, the transport of solutes remains limited.

4.2. Simulations of Individual Glacial Processes

4.2.1. Simulations B1 and B2: Subglacial Recharge Boundary Condition

In simulation B1, where a 2 mm/y flux is assigned, the constraint does not activate because the hydraulic heads never reach the pressure of the ice-sheet. Meltwater transport is limited to the upper portions of the aquifers. Groundwater flow velocities do not exceed 250 mm/y and meltwater content reaches 5 and 4% in total in the intermediate and lower aquifers.

In simulation B2, where a 10 mm/y flux is assigned, meltwater inflow is sufficient for hydraulic head to quickly reach the pressure of the ice-sheet. Accordingly, the prescribed flux boundary condition turns into a fixed head boundary condition. Groundwater flow and the transport of solutes evolve in a similar way to reference simulation A.

4.2.2. Simulations C: Glacial Loading

Direct loading equilibrates instantaneously the hydraulic head with surface conditions. It significantly decreases abnormal pressures in confining layers: in the lower aquitard, overpressures are about 270 and 40 m after the final retreat of the ice-sheet and at the end of the simulation, against 800 and 70 m in reference simulation A. Instantaneous change of hydraulic head within the lower aquitard decreases the flow rates in the intermediate and lower aquifers, which is illustrated by less meltwater infiltration. In the rest of the basin, hydraulic diffusivity is too high for direct loading to have any impact on hydraulic head.

4.2.3. Simulations D1 and D2: Isostasy

As showed on Figure 7, the northern part of the basin is lowered with respect to the southern part, due to isostasy. It results in a topographic gradient opposite to the ice-sheet gradient. In consequence, the resulting hydraulic gradient beneath the ice-sheet is smaller and less meltwater infiltrates compared to reference simulation A. During the ice-sheet retreat, the topographic gradient remains in simulation D1 due to the postglacial isostatic depression, which decreases from 80 to 3 m in the northern part of the basin. In simulation D2, this topographic gradient is cancelled by the water body that fills up the postglacial depression. Accordingly, more meltwater is expected to flow out of the basin during the postglacial period in simulation D1 with respect to simulation D2. However, this effect is barely noticeable.

4.2.4. Simulation E: Permafrost

Figure 8 shows the degradation of permafrost during the growth of the ice-sheet. The presence of permafrost delays the infiltration of meltwater into the subsurface of about 1 ky. As a consequence, the maximal meltwater content reached in each aquifer is lower than in the reference simulation. Flow rates below the permafrost layer are lower than in reference simulation A but they are enhanced further beyond the front of the ice-sheet.

4.2.5. Simulation F: Variable Density

In the lower part of the basin, hydraulic heads are naturally larger, because denser water increases the hydrostatic pressure. Accordingly, the subglacial regional gradient is less than in reference simulation A, especially in the lowest aquifer. Deeper than the outcrop zone of the lowest aquifer, flow rates are about 50% lower than in reference simulation A. In consequence, less meltwater infiltrates. Meltwater content is about 15% less than in the reference simulation.

5. Discussion

In this model, representing subglacial recharge by a fixed head or a fixed flux constrained by the same fixed head has no influence on the solution for values of fixed flux equal to or larger than 10 mm/y. Such flux guarantees a fast increase of hydraulic heads and activation of the constraint. The time necessary for hydraulic heads to reach the constraint depends on the hydraulic diffusivity of the subsurface and on the value of the constraint (i.e., the head equivalent to the weight of the ice-sheet). In the present model, hydraulic diffusivity and the upper constraint have intermediate values. It follows that a fixed head boundary condition is a suitable representation of subglacial recharge for many natural cases. Fixed flux should be used only if the flux is known to be small (not larger than a few mm/y), which could be the case for partially wet-based ice-sheets, or if the hydraulic diffusivity and the upper constraint have, respectively, small and large values.

Direct loading considerably reduces abnormal pressures in the basin. Therefore, it has an impact only if abnormal pressures are generated beneath the ice-sheet. In this study, overpressures are noticed only in the lower aquitard, where the hydraulic diffusivity is about 10⁻⁶ m²/s. This range of values is really low but consistent with the formula estimating the time required for dissipating abnormal pressures in hydrogeological units:where is half the thickness of the layer [m] [33]. As in Black and Barker [50], this model reproduces overpressures simulating subglacial recharge, with or without direct loading. To simulate underpressures in the Michigan Basin, Nasir et al. [34], Khader and Novakowski [25], and Neuzil and Provost [35] assigned values of subglacial recharge boundary conditions inferior to the value of ice-sheet loading, in order to allow overpressures induced by glacial loading to drain out at the surface and then turn into underpressures. Some of these models even assumed no subglacial recharge at all, which would correspond to a dry-based ice-sheet scenario. It follows that abnormal pressures are dependent not only on direct loading but also on subglacial recharge, which was already acknowledged by Khader and Novakowski [25].

During the advance of the ice-sheet, subglacial gradient and flow rates are reduced by a percentage similar to the ratio between the ice-sheet thickness and the isostatic depletion. Therefore, simulating subglacial recharge beneath ice-sheets without considering isostasy could lead to a major overestimation of subglacial flow. More meltwater was expected to be preserved in the simulation considering postglacial lake or sea. Here, it is likely that this effect is masked by the dissipation of overpressures from the lower aquitard, which has a larger impact on groundwater flow. Studies showed that large water bodies preserve meltwater from postglacial topography-driven recharge [12, 21, 26, 43]. Given the flat topography assigned in the present model, this effect also passes unnoticed.

As noted in the studies cited in Section 2.4, permafrost does prevent groundwater from being flushed out through the surface. In consequence, high groundwater pressures propagate beyond the front of the ice-sheet and groundwater flow rates are smoothed near the front of the ice-sheet. In the present study, a smooth profile of hydraulic heads was assigned beneath the ice-sheet so that permafrost does not influence the maximal flow rate. If a peak of groundwater flow rates happens at the front of the ice-sheet, as in Vidstrand et al. [39], permafrost could have such an influence. Permafrost also has a retardation effect on the infiltration of meltwater into the subsurface. An extreme case would be the nonthawing of permafrost beneath a dry-based ice-sheet, as simulated in Nasir et al. [34] and Khader and Novakowski [25]. In this case, subglacial recharge would be limited to taliks, if present, or glaciations would correspond to gaps in recharge [51].

Variable density simulation shows that the flow rates and the infiltration of meltwater are restrained by larger hydrostatic heads at depth. However, large volumes of meltwater still infiltrate, which is in agreement with previous models that also considered the variations in density [18, 19, 22, 24, 27]. The presence of heavy brines slows down the displacement of meltwater into the deeper part of the basin but does not restrict it to the edge of the basin.

This study shows that every glacial process has a potentially large influence on groundwater flow in terms of meltwater infiltration and groundwater flow rates. The only process that could be discarded a priori is direct loading, if thick and low permeability aquitards are absent. However, some glacial processes could be neglected in order to produce models that are easier to run. Given the uncertainty related to the representation of glacial processes under past and a fortiori future environmental conditions, many variables arise. Then, neglecting some glacial processes could also serve to avoid problems related to nonuniqueness of simulation results. The results of this study provide a good basis to decide which glacial processes to neglect, if necessary, and to determine the impact of neglecting them. However, the scope of this study has some limits. In natural basins, different geometries could lead to different results when simulating glacial processes. There is also a possibility that simulating different glacial processes together could be nonlinear; that is, the respective effects of each process could not be simply summed up. In consequence, this study should not replace a proper sensitivity analysis in future hydrogeological modelling studies of natural systems.

6. Conclusion

This study assessed the specific impact of glacial processes associated with wet-based ice-sheets that were simulated in previous hydrogeological modelling studies: subglacial recharge, direct glacial loading, isostasy, the presence of proglacial lake or sea, permafrost, and density-dependent mixing between meltwater and brines. Those processes were incorporated in a groundwater flow and solute transport numerical model of a generic sedimentary basin. Several simulations were performed in order to assess the impact of these glacial processes in terms of meltwater infiltration and groundwater flow rates. Results show that all processes are potentially relevant with respect to groundwater flow and clarify how relevant they are.

Assigning a fixed flux boundary condition constrained by a fixed head equivalent to the weight of the ice-sheet is meaningful only if the flux of meltwater is low and does not allow hydraulic head to reach the constraint. In the present model, where intermediate values of hydraulic diffusivity and constraint were applied, fluxes equal to or larger than 10 mm/y show no difference with fixed head boundary conditions.

Direct glacial loading has an impact only if abnormal pressures are generated, which happens only in thick and low permeability aquitards. Then, direct loading decreases overpressures created by subglacial recharge. If subglacial recharge is less than direct loading, underpressures appear. Equation (9) can be used to estimate which hydrogeological units will host abnormal pressures, on the basis of their thickness and their hydraulic diffusivity.

Isostasy lowers the part of the basin under the ice-sheet with respect to the rest of the basin. It creates a hydraulic gradient opposite to the direction of the ice-sheet. In consequence, it acts against the infiltration of meltwater and reduces the groundwater flow rates. The preservation of meltwater under postglacial water bodies passes unnoticed because of the dissipation of overpressures from the lower aquitard and the flat topography used in the model.

Permafrost prevents the subsurface to be recharged by subglacial meltwater. If the thawing time is long, the volume of meltwater that is recharged during a glaciation can be considerably reduced. Beneath permafrost, hydraulic heads increase beyond the ice-sheet front. Accordingly, groundwater flow is reduced in shallow aquifers beneath the ice-sheet but is enhanced beyond the ice-sheet front.

Accounting for salinity-dependent density increases hydrostatic pressure in the deeper part of the basin. The infiltration of meltwater at depth is reduced.

In conclusion, none of these glacial processes can be discarded a priori, except direct loading in cases. However, some could be neglected purposely to produce simpler models, especially if information is missing to describe them properly. In this case, this study provides useful means to determine the impact of such simplification.

Disclosure

Earlier versions of this work were presented in an oral presentation at the G@GPS Project International Workshop in Tallinn (2015) and as a poster at the American Geophysical Union Fall Meeting in San Francisco (2015).

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was supported by the Discovery Grant from the Natural Sciences and Engineering Research Council of Canada attributed to Jean-Michel Lemieux, by the Estonian Research Council Project IUT 19-22 and by the DoRa Programme.

References

D. I. Siegel and R. J. Mandle, “Isotopic evidence for glacial meltwater recharge to the Cambrian-Ordovician aquifer, north-central United States,” Quaternary Research, vol. 22, no. 3, pp. 328–335, 1984.
View at: Publisher Site | Google Scholar
I. D. Clark, M. Douglas, K. Raven, and D. Bottomley, “Recharge and preservation of Laurentide glacial melt water in the Canadian shield,” Ground Water, vol. 38, no. 5, pp. 735–742, 2000.
View at: Publisher Site | Google Scholar
R. Vaikmäe, L. Vallner, H. H. Loosli, P. C. Blaser, and M. Juillard-Tardent, “Palaeogroundwater of glacial origin in the Cambrian-Vendian aquifer of Northern Estonia,” Geological Society Special Publication, vol. 189, pp. 17–27, 2001.
View at: Publisher Site | Google Scholar
J. C. Mcintosh, M. E. Schlegel, and M. Person, “Glacial impacts on hydrologic processes in sedimentary basins: evidence from natural tracer studies,” Geofluids, vol. 12, no. 1, pp. 7–21, 2012.
View at: Publisher Site | Google Scholar
G. S. Boulton, P. E. Caban, and K. Van Gijssel, “Groundwater flow beneath ice sheets: part I—large scale patterns,” Quaternary Science Reviews, vol. 14, no. 6, pp. 545–562, 1995.
View at: Publisher Site | Google Scholar
J. Kleman and N. F. Glasser, “The subglacial thermal organisation (STO) of ice sheets,” Quaternary Science Reviews, vol. 26, no. 5-6, pp. 585–597, 2007.
View at: Publisher Site | Google Scholar
G. S. Boulton, P. E. Caban, K. Van Gijssel, A. Leijnse, M. Punkari, and F. H. A. Van Weert, “The impact of glaciation on the groundwater regime of Northwest Europe,” Global and Planetary Change, vol. 12, no. 1–4, pp. 397–413, 1996.
View at: Publisher Site | Google Scholar
S. E. Grasby and Z. Chen, “Subglacial recharge into the Western Canada Sedimentary Basin—impact of Pleistocene glaciation on basin hydrodynamics,” Bulletin of the Geological Society of America, vol. 117, no. 3-4, pp. 500–514, 2005.
View at: Publisher Site | Google Scholar
T. Grundl, N. Magnusson, M. S. Brennwald, and R. Kipfer, “Mechanisms of subglacial groundwater recharge as derived from noble gas, 14C, and stable isotopic data,” Earth and Planetary Science Letters, vol. 369-370, pp. 78–85, 2013.
View at: Publisher Site | Google Scholar
S. Grasby, K. Osadetz, R. Betcher, and F. Render, “Reversal of the regional-scale flow system of the Williston basin in response to Pleistocene glaciation,” Geology, vol. 28, no. 7, pp. 635–638, 2000.
View at: Publisher Site | Google Scholar
C. W. Breemer, P. U. Clark, and R. Haggerty, “Modeling the subglacial hydrology of the late Pleistocene Lake Michigan Lobe, Laurentide Ice Sheet,” Bulletin of the Geological Society of America, vol. 114, no. 6, pp. 665–674, 2002.
View at: Publisher Site | Google Scholar
M. Person, B. Dugan, J. B. Swenson et al., “Pleistocene hydrogeology of the Atlantic continental shelf, New England,” Bulletin of the Geological Society of America, vol. 115, no. 11, pp. 1324–1343, 2003.
View at: Publisher Site | Google Scholar
C. Talbot, “Ice ages and nuclear waste isolation,” Engineering Geology, vol. 52, no. 3-4, pp. 177–192, 1999.
View at: Publisher Site | Google Scholar
M. F. Loutre and A. Berger, “Future climatic changes: are we entering an exceptionally long interglacial?” Climatic Change, vol. 46, no. 1-2, pp. 61–90, 2000.
View at: Publisher Site | Google Scholar
M. Person, V. Bense, D. Cohen, and A. Banerjee, “Models of ice-sheet hydrogeologic interactions: a review,” Geofluids, vol. 12, no. 1, pp. 58–78, 2012.
View at: Publisher Site | Google Scholar
I. Lerche, Z. Yu, B. Tørudbakken, and R. O. Thomsen, “Ice loading effects in sedimentary basins with reference to the Barents Sea,” Marine and Petroleum Geology, vol. 14, no. 3, pp. 277–338, 1997.
View at: Publisher Site | Google Scholar
J. C. McIntosh and L. M. Walter, “Volumetrically significant recharge of Pleistocene glacial meltwaters into epicratonic basins: constraints imposed by solute mass balances,” Chemical Geology, vol. 222, no. 3-4, pp. 292–309, 2005.
View at: Publisher Site | Google Scholar
J.-M. Lemieux, E. A. Sudicky, W. R. Peltier, and L. Tarasov, “Simulating the impact of glaciations on continental groundwater flow systems: 2. Model application to the Wisconsinian glaciation over the Canadian landscape,” Journal of Geophysical Research: Earth Surface, vol. 113, no. 3, Article ID F03018, 2008.
View at: Publisher Site | Google Scholar
V. F. Bense and M. A. Person, “Transient hydrodynamics within intercratonic sedimentary basins during glacial cycles,” Journal of Geophysical Research: Earth Surface, vol. 113, no. 4, 2008.
View at: Publisher Site | Google Scholar
J.-M. Lemieux, E. A. Sudicky, W. R. Peltier, and L. Tarasov, “Simulating the impact of glaciations on continental groundwater flow systems: 1. Relevant processes and model formulation,” Journal of Geophysical Research: Earth Surface, vol. 113, no. 3, Article ID F03017, 2008.
View at: Publisher Site | Google Scholar
D. Cohen, M. Person, P. Wang et al., “Origin and extent of fresh paleowaters on the Atlantic continental shelf, USA,” Ground Water, vol. 48, no. 1, pp. 143–158, 2010.
View at: Publisher Site | Google Scholar
J. C. McIntosh, G. Garven, and J. S. Hanor, “Impacts of Pleistocene glaciation on large-scale groundwater flow and salinity in the Michigan Basin,” Geofluids, vol. 11, no. 1, pp. 18–33, 2011.
View at: Publisher Site | Google Scholar
S. D. Normani and J. F. Sykes, “Paleohydrogeologic simulations of Laurentide ice-sheet history on groundwater at the eastern flank of the Michigan Basin: paleohydrogeologic simulations on the eastern flank of the Michigan Basin,” Geofluids, vol. 12, no. 1, pp. 97–122, 2012.
View at: Publisher Site | Google Scholar
A. M. Provost, C. I. Voss, and C. E. Neuzil, “Glaciation and regional groundwater flow in the Fennoscandian shield,” Geofluids, vol. 12, no. 1, pp. 79–96, 2012.
View at: Publisher Site | Google Scholar
O. Khader and K. Novakowski, “Impacts of Pleistocene glacial loading on abnormal pore-water pressure in the eastern Michigan Basin,” Geofluids, vol. 14, no. 2, pp. 200–220, 2014.
View at: Publisher Site | Google Scholar
J. Siegel, M. Person, B. Dugan, D. Cohen, D. Lizarralde, and C. Gable, “Influence of late Pleistocene glaciations on the hydrogeology of the continental shelf offshore Massachusetts, USA,” Geochemistry, Geophysics, Geosystems, vol. 15, no. 12, pp. 4651–4670, 2014.
View at: Publisher Site | Google Scholar
P. Vidstrand, S. Follin, J.-O. Selroos, and J.-O. Näslund, “Groundwater flow modeling of periods with periglacial and glacial climate conditions for the safety assessment of the proposed high-level nuclear waste repository site at Forsmark, Sweden,” Hydrogeology Journal, vol. 22, no. 6, pp. 1251–1267, 2014.
View at: Publisher Site | Google Scholar
J. A. Piotrowski, “Subglacial hydrology in north-western Germany during the last glaciation: groundwater flow, tunnel valleys and hydrological cycles,” Quaternary Science Reviews, vol. 16, no. 2, pp. 169–185, 1997.
View at: Publisher Site | Google Scholar
H. Engelhardt and B. Kamb, “Basal hydraulic system of a West Antarctic ice stream: constraints from borehole observations,” Journal of Glaciology, vol. 43, no. 144, pp. 207–230, 1997.
View at: Publisher Site | Google Scholar
L. Claesson Liljedahl, A. Kontula, J. Harper et al., “The Greenland Analogue Project: final report,” Tech. Rep. SKB TR-14-13, Svensk Kärnbränslehantering, Stockholm, Sweden, 2016.
View at: Google Scholar
V. J. Lunardini, “Climatic warming and the degradation of warm permafrost,” Permafrost and Periglacial Processes, vol. 7, no. 4, pp. 311–320, 1996.
View at: Publisher Site | Google Scholar
F. H. A. van Weert, K. van Gijssel, A. Leijnse, and G. S. Boulton, “The effects of Pleistocene glaciations on the geohydrological system of Northwest Europe,” Journal of Hydrology, vol. 195, no. 1–4, pp. 137–159, 1997.
View at: Publisher Site | Google Scholar
C. E. Neuzil, “Hydromechanical effects of continental glaciation on groundwater systems,” Geofluids, vol. 12, no. 1, pp. 22–37, 2012.
View at: Publisher Site | Google Scholar
O. Nasir, M. Fall, T. S. Nguyen, and E. Evgin, “Modelling of the hydro-mechanical response of sedimentary rocks of southern Ontario to past glaciations,” Engineering Geology, vol. 123, no. 4, pp. 271–287, 2011.
View at: Publisher Site | Google Scholar
C. E. Neuzil and A. M. Provost, “Ice sheet load cycling and fluid underpressures in the Eastern Michigan Basin, Ontario, Canada,” Journal of Geophysical Research B: Solid Earth, vol. 119, no. 12, pp. 8748–8769, 2014.
View at: Publisher Site | Google Scholar
W. Rühaak, V. F. Bense, and I. Sass, “3D hydro-mechanically coupled groundwater flow modelling of Pleistocene glaciation effects,” Computers & Geosciences, vol. 67, pp. 89–99, 2014.
View at: Publisher Site | Google Scholar
R. L. Kleinberg and D. D. Griffin, “NMR measurements of permafrost: unfrozen water assay, pore-scale distribution of ice, and hydraulic permeability of sediments,” Cold Regions Science and Technology, vol. 42, no. 1, pp. 63–77, 2005.
View at: Publisher Site | Google Scholar
K. Horiguchi and R. D. Miller, “Experimental studies with frozen soil in an “ice sandwich” permeameter,” Cold Regions Science and Technology, vol. 3, no. 2-3, pp. 177–183, 1980.
View at: Publisher Site | Google Scholar
P. Vidstrand, S. Follin, J.-O. Selroos, J.-O. Näslund, and I. Rhén, “Modeling of groundwater flow at depth in crystalline rock beneath a moving ice-sheet margin, exemplified by the Fennoscandian Shield, Sweden,” Hydrogeology Journal, vol. 21, no. 1, pp. 239–255, 2013.
View at: Publisher Site | Google Scholar
J. Harff and M. Meyer, “Coastlines of the Baltic Sea—zones of competition between geological processes and a changing climate: examples from the Southern Baltic,” in The Baltic Sea Basin, pp. 149–164, Springer, Berlin, Germany, 2011.
View at: Publisher Site | Google Scholar
S. Björck, “A review of the history of the Baltic Sea, 13.0–8.0 ka BP,” Quaternary International, vol. 27, pp. 19–40, 1995.
View at: Publisher Site | Google Scholar
G. K. C. Clarke, D. W. Leverington, J. T. Teller, and A. S. Dyke, “Paleohydraulics of the last outburst flood from glacial Lake Agassiz and the 8200 BP cold event,” Quaternary Science Reviews, vol. 23, no. 3-4, pp. 389–407, 2004.
View at: Publisher Site | Google Scholar
W. DeFoor, M. Person, H. C. Larsen, D. Lizarralde, D. Cohen, and B. Dugan, “Ice sheet-derived submarine groundwater discharge on Greenland's continental shelf,” Water Resources Research, vol. 47, no. 7, 2011.
View at: Publisher Site | Google Scholar
S. A. Bea, U. K. Mayer, and K. T. B. MacQuarrie, “Reactive transport and thermo-hydro-mechanical coupling in deep sedimentary basins affected by glaciation cycles: model development, verification, and illustrative example,” Geofluids, vol. 16, no. 2, pp. 279–300, 2016.
View at: Publisher Site | Google Scholar
C. van der Veen, Fundamentals of Glacier Dynamics, CRC Press, Boca Raton, Fla, USA, 2013.
View at: Publisher Site
R. L. Stotler, S. K. Frape, T. Ruskeeniemi, P. Pitkänen, and D. W. Blowes, “The interglacial-glacial cycle and geochemical evolution of Canadian and Fennoscandian Shield groundwaters,” Geochimica et Cosmochimica Acta, vol. 76, pp. 45–67, 2012.
View at: Publisher Site | Google Scholar
P. Vidstrand, U. Svensson, and S. Follin, “Simulation of hydrodynamic effects of salt rejection due to permafrost,” Tech. Rep. SKB R-06-101, Svensk Kärnbränslehantering, Stockholm, Sweden, 2006.
View at: Google Scholar
J. S. Hanor, “Sedimentary genesis of hydrothermal fluids,” in Geochemistry of Hydrothermal Ore Deposits, H. L. Barnes, Ed., pp. 137–168, John Wiley & Sons, New York, NY, USA, 1979.
View at: Google Scholar
H.-J. G. Diersch, FEFLOW, Springer, Berlin, Germany, 2014.
J. H. Black and J. A. Barker, “The puzzle of high heads beneath the West Cumbrian coast, UK: a possible solution,” Hydrogeology Journal, vol. 24, no. 2, pp. 439–457, 2016.
View at: Publisher Site | Google Scholar
H. Jiráková, F. Huneau, H. Celle-Jeanton, Z. Hrkal, and P. Le Coustumer, “Insights into palaeorecharge conditions for European deep aquifers,” Hydrogeology Journal, vol. 19, no. 8, pp. 1545–1562, 2011.
View at: Publisher Site | Google Scholar

Copyright

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

PDF Download Citation

Download other formats

Order printed copies

Views

1689

Downloads

1419

Citations