Changing groundwater discharge dynamics in permafrost regions

We use a modified version of the USGS finite element groundwater model SUTRA, which simulates groundwater flow in variably-saturated conditions with conductive-advective energy transport (Voss and Provost 2010). As described in McKenzie et al (2007), the present version of SUTRA includes pore-water phase change and the resultant impacts on permeability, fluid density and release or absorption of heat due to the latent heat of fusion. It has been successfully tested against numerical and analytical solutions (Kurylyk et al 2014b), including the international benchmarking initiative, InterFrost (Rühaak et al 2015, Grenier et al 2018), and used in many permafrost studies to model groundwater flow dynamics (e.g. Ge et al 2011, McKenzie and Voss 2013, Briggs et al 2014, Kurylyk et al 2016, Evans and Ge 2017, Zipper et al 2018). The present code version allows for unsaturated flow with freeze/thaw as described in Briggs et al (2014).

The intrinsic permeability of the unsaturated and/or frozen ground is a function of liquid water content. When saturated and unfrozen, the horizontal permeability is 10⁻¹³ m² and exponentially reduces to 10⁻¹⁹ m² at the minimum liquid water content. These values are representative of permafrost environments and have been used in previous modelling studies (Bense et al 2012, McKenzie and Voss 2013, Evans and Ge 2017). Pore water begins freezing when the temperature is below 0 °C, and the pore-ice saturation increases with decreasing temperature to −2 °C via an exponential soil freezing curve (Mottaghy and Rath 2006). See the supporting information available at stacks.iop.org/ERL/13/084017/mmedia for further information.

The model domain is a 500 m long, two-dimensional homogeneous cross-section with a 2% land surface slope (figure 1(a)). This model shape represents a simple continuous permafrost landscape where groundwater discharges into a medium sized river. One half of a triangular river channel cross-section is represented by a 50% land-surface slope for the last 10 m of the cross section (inset, figure 1(a)). The model thickness is 410 m upslope and 395 m downslope. Variations of this domain geometry were considered in a sensitivity study as described later. The vertical node spacing increases from 0.5 m at the land surface to 10 m at the model bottom, and the horizontal node spacing is 5.0 m. To ensure that the grid cell and time step size did not impact our results, we performed a mesh refinement sensitivity analysis by comparing the groundwater discharge results and the total volume of ice in the model. Smaller time steps or node spacings influenced the quantity of ice in the model by less than 0.04%, giving us confidence that the discretization did not impact our results. The mesh refinement study (results in the supplementary documentation) is different from the sensitivity analysis as it only measures the effects of changing the model resolution.

The vertical boundaries of the model represent hydrogeological divides with no heat flux or water flow. The bottom edge of the domain is a no-flow hydrologic boundary since model runs indicated that there is no significant groundwater flow below 400 meters depth. A geothermal heat flux of 0.085 W m⁻² is assigned across the bottom of the domain (McKenzie and Voss 2013). The domain upper surface has two sections: the rightmost 10 m represents a river, and the remainder of the top boundary represents the land surface (figure 1(a)). The land surface is represented by a superimposed combination of specified recharge, hydrologic drain, and specified heat flux BCs. Using a specified recharge instead of a specified pressure BC, as in most previous studies (e.g. Ge et al 2011, Bense et al 2012, McKenzie and Voss 2013), better represents groundwater recharge in a climate-controlled system (Sanford 2002).

The specified groundwater recharge boundary is a function of the time of the year and climatic conditions (i.e. air temperature (AT), summer precipitation (SP) and winter precipitation (WP)). The AT is defined as a sinusoidal function with a seasonal temperature amplitude (i.e. half the total range) of 19.5 °C (figure 1(b)). The SP and WP are 56.5 mm month⁻¹ and 24.3 mm month⁻¹, respectively, based on averages from Churchill, Canada which has continuous permafrost (58°46'N 094°10'W; Environment Canada 2017). The seasonal specified recharge BC has 3 stages (figure 1(b)). First, when AT is below 0 °C, recharge is zero and WP accumulates as snowpack on the land surface. The second stage simulates freshet—for the 15 days after AT first exceeds 0 °C, recharge rapidly increases, representing the spring freshet. Recharge during this period comes from both SP and snowmelt, and the peak value is proportional to the amount of snow accumulated during the below freezing season (BFS, figure 1(b)). During the 15 days following this peak, recharge decreases linearly until no snowpack remains. Third, after freshet, recharge is equal to 20% of SP. Once AT falls below 0 °C, the annual cycle repeats. This BC is representative of groundwater recharge rates observed in field measurements in permafrost environments (e.g. Clilverd et al 2011, Murray 2016). This BC also allow for change in annual precipitation: as mean annual AT increases, the duration of summer increases and, since SP is higher than WP, the model receives an increase of annual precipitation.

The surface drain BC allows over-pressured surface nodes to discharge water. The quantity of water being discharged from the drain nodes is a function of the simulated pressure at the node. Unlike a specified pressure boundary, this combination of specified recharge and drain conditions prevents an unrealistic amount of recharge from being forced into the model at the land surface. Furthermore, these boundaries allow groundwater to exit the model anywhere along the land surface rather than forcing discharge to only occur at specified locations (e.g. a river boundary) as most prior studies have assumed. Because it is applied in conjunction with a drain, the specified recharge BC defines the maximum available water for recharge, not the actual recharge to the aquifer.

Zoom In Zoom Out Reset image size

For the land surface thermal BC, previous modelling studies typically used either a specified ground surface temperature (GST) boundary (Bense et al 2012, Frampton et al 2013) or a specified AT with an underlying thermal boundary layer (McKenzie et al 2007, Evans and Ge 2017). The former may negate the effect of subsurface processes on the surface temperatures by prescribing GST, while the latter ignores the insulating effect of winter snowpack. To overcome these limitations, we use a specified heat flux boundary at the land surface calculated with a conduction equation based on the temperature difference between the specified AT and the modelled GST. This novel BC accounts for changes in the thermal conductivity between the air and the surface due to the insulating effect of snow cover (supporting information).

The riverbed is simulated using a specified pressure boundary representing a river with the pressure at the river-bottom set as hydrostatic in accordance with the water depth. The thermal boundary at the river is a variable specified temperature boundary that mimics the temperatures at the bottom of a wide river or a lake. The temperatures are calculated from a synthetic function based on lake-bottom temperatures measured in a permafrost environment (Burn 2003) and depend on the time of the year, the AT and river-bottom depth. The function generates results similar to that of Wellman et al (2013). This simplified river representation is used to isolate groundwater flow processes which are the focus of the present study.

The model simulation is run in three phases; phases 1 and 2 generate the initial temperature, permafrost, and pressure distributions that are used as the starting point for phase 3 that produces the simulation results. Phase 1 generates a permafrost and temperature distribution by simulating 30 000 years with a constant air temperature of −6.5 °C without seasonal variation or recharge at the land surface. At the end of the phase 1, the model is at steady state with a mean permafrost depth of ∼250 m. Phase 2 creates a quasi-dynamic equilibrium for the near-surface permafrost and active zone by running the model for an additional 100 years with seasonal recharge and an annually varying surface temperature, with a mean of −6.5 °C and an amplitude of 19.5 °C (figure 1(b)). During phase 3, AT increases in a piecewise linear fashion, superimposed with a sinusoidal annual air temperature fluctuation (figure 1(c)). From 0 to 100 years, the warming rate is 0.06 °C year⁻¹ and represents the average rate of increase from multi-model climate projections for the Arctic (Kattsov et al 2005). From year 100–300, the warming rate is 0.02 °C year⁻¹, which is a mix between short-term projections for the Arctic and global long-term projections (Kattsov et al 2005, Zickfeld et al 2013). After 300 years, the warming rate is 0.0015 °C year⁻¹ and represents the average of long-term global climate projections (Zickfeld et al 2013). The time step sizes are 10 years for phase 1 and 2 hours for phases 2 and 3. For the results presented below, time = 0 corresponds to the beginning of phase three.

In order to demonstrate that the groundwater discharge variations and patterns simulated in this study are not limited to one specific model configuration, eight additional simulations are used. Permeability, surface slope, recharge rate and snow density (winter insulation) have been selected as parameters that may have an important impact on groundwater discharge and permafrost distribution (e.g. McKenzie and Voss 2013). For each simulation, one parameter is given either a lower or a higher value than the one specified in the default simulation, while the other parameters remain constant (table 1).

Figure 2 conceptually illustrates the evolution of permafrost distribution and groundwater flow through different times of the climate warming simulation. A 14 m deep and 12 m wide talik is initially present under the simulated river (figure 2(a)). During the first 70 years of warming, a mushy zone (i.e. where both liquid water and ice coexist, Amiri et al 2018) gradually forms at the interface between the active layer and permafrost (figures 2(b), 3(a)), which increases the permeability and enables slow groundwater flow even below the permafrost table (figure 3(b)). The mushy zone in figure 3(a) shows the area that stays above −0.25 °C all year, meaning that the permeability never decreases more than one order of magnitude. This mushy zone opens to create a permanently unfrozen zone (talik) between the active layer and the permafrost table after approximately 70 years (figures 2(c), 3(a)). With warming, the summers become longer and warmer, leading to deeper annual thawing. The resultant shorter winters are not sufficiently cold and/or long to completely counteract this thawing. As the horizontal talik extends, the increased groundwater flow (figure 3(b)) has a positive feedback on the talik's growth as year-round heat advection becomes more important and drives further thawing. Heat advection in the talik during the winter season, combined with warmer winter air temperatures, causes a decrease of the seasonal frost depth, and thus the active layer thickness, once the talik starts forming after 70 years (figure 3(a)). After 200 years of warming, the supra-permafrost aquifer reaches a sufficient thickness to accommodate the total annual recharge. This increased groundwater storage is evidenced by the fact that the water table remains primarily below the land surface throughout the entire year (figure 2(d)). For the given climate warming rates, all permafrost thaws after 992 years.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Groundwater discharge can occur via two pathways: (1) land surface drain nodes, subdivided into upgradient (0 < × < 245 m) and downgradient (245 < × < 490) zones, and (2) specified pressure nodes representing the riverbed (figure 1(a)). There is a pronounced shift in the spatial distribution of groundwater discharge as permafrost thaws (figure 4(a)). Initially, most of the groundwater discharges to the land surface as opposed to the river because the pore ice limits connections between the river and aquifer. This discharge is approximately evenly distributed between upgradient and downgradient areas. Such groundwater discharge settings are typical of poorly drained wetlands systems in permafrost environments (Woo 2012). As the system warms, groundwater discharge to the upgradient land surface decreases drastically, and discharge shifts to the river boundary. After 300 years, the annual volume of water discharging through the upgradient surface decreases by 50%, and the annual volume of water discharging at the riverbed is 15 times higher than at the beginning (figure 4(a)).

Zoom In Zoom Out Reset image size

To quantify the groundwater discharge distribution shift, we use the GWDDR, calculated as the ratio of specific annual groundwater discharge through the riverbed over specific annual groundwater discharge through the land surface. Increasing supra-permafrost thickness is the main factor driving the increase in GWDDR (figure 4(b), shift in discharge to the riverbed) until the aquifer becomes thick enough to accommodate the total recharge, after which GWDDR remains relatively constant (figure 4(a)).

We suggest that these shifts in the spatial distribution of groundwater discharge may explain the paradox of why some wetlands and lakes have been observed to drain (e.g. Smith et al 2005, Necsoiu et al 2013), yet baseflow to rivers has increased in permafrost environments (e.g. Walvoord and Striegl 2007, Déry et al 2009, St Jacques and Sauchyn 2009). Our heuristic modeling approach allows us to identify the processes underlying this spatial shift in groundwater discharge. As permafrost thaws, new lateral pathways form (figure 2(d)), which increase hydrogeological connectivity between the recharge zones and the riverbed. Conversely when shallow permafrost is present, it forces water to discharge and exfiltrate the hillslope aquifer via the drain boundary along the land surface (figure 2(a)). It should be noted that these variations occur during the thawing of only 40% of the total permafrost, suggesting that shallow permafrost thaw can have a disproportionately large effect on the spatial distribution of groundwater discharge and surface water features. This provides mechanistic support for Jepsen et al (2013), who proposed that changes in lake area are driven by changes in shallow groundwater flow paths. Furthermore, the results indicate that field or remote sensing observations of disappearing wetlands in permafrost terrains may indicate lateral talik formation.

To analyse the seasonality of groundwater discharge, the results are separated into two periods: the below freezing season (BFS), when the AT is below 0 ˚C, and the above freezing season (AFS), when the AT is above 0 ˚C. As mean AT increases, the BFS shortens and the AFS lengthens. Because SP from the specified recharge BC is higher than WP, the model receives a yearly increase in groundwater recharge during the warming phase. After 100 years, the cumulative increase is 9%, which is within the range of arctic precipitation projections (Kattsov et al 2005).

Zoom In Zoom Out Reset image size

Figures 5(a) and (b) compares the mean daily outflows (i.e. discharge normalized to boundary length) during these periods. The BFS discharge increases throughout the entire warming simulation, with the largest increase between 40 and 100 years. This early rapid increase is explained by the formation of a talik between the bottom of the frozen active layer and the permafrost table (figure 2(c)), which provides a winter groundwater connection between upslope areas and the unfrozen riverbed. The BFS discharge occurs almost entirely through the riverbed (figure 5(b)), which highlights the impact of the increasing connectivity due to the talik formation on discharge variations.

Although field studies about lateral supra-permafrost taliks are scarce, recent publications have reported them in both discontinuous (Zhou et al 2015, Connon et al 2018) and continuous (Hiyama et al 2013) permafrost regions and primarily focused on the climate and environmental controls of their occurrence (e.g. Connon et al 2018). While Hiyama et al (2013) provided some observations of talik derived groundwater discharge, no field studies have yet quantified the importance of lateral taliks to groundwater-surface water interactions. Given that our modeling results indicate lateral taliks can be critical controls on surface-subsurface hydrologic connectivity, additional field studies are needed to determine where and under what conditions these features are important for arctic hydrology.

The AFS discharge evolution exhibits the opposite pattern (figure 5(a)). The downward migration of the permafrost table enables the water table to stabilize, reducing the hydraulic gradient of the aquifer and the concomitant summer groundwater discharge (figure 2). However, there is no significant change in the water table after 200 years indicating that factors besides the decrease in hydraulic gradient are important. Rather, the permafrost thaw and lateral talik formation allows more liquid water to be stored in the groundwater system throughout the year (figure 2), thus reducing the groundwater discharge seasonality.

Zoom In Zoom Out Reset image size

We use the GWDSR, defined as the ratio of the mean annual groundwater discharge during AFS over the mean annual groundwater discharge during BFS to quantify temporal shifts in groundwater discharge. Figure 5(c) shows that the decrease in seasonality is initially caused by a longer duration of seasonally thawed ground. Once this zone becomes permanent (creating a lateral talik), increasing talik thickness leads to a decrease in groundwater discharge seasonality.

Our results suggest that the measured increase of arctic river baseflow may be due to a combination of increasing precipitation and connectivity in aquifers. Since the groundwater discharge remains relatively constant after only 40% of the total permafrost thaws, the contribution of meltwater from permafrost is not an important factor causing increased baseflow. Like discharge, recharge rates approach a steady rate after 300 years. However, the simulation results show that the early variations are too high to be explained by an increase in precipitation alone (which is 9% after 100 years in this simulation). Rather, baseflow changes are predominantly driven by the hydrogeological connectivity increase due to shallow permafrost thaw.

While our results demonstrate the importance of increased subsurface connectivity, other factors may also contribute to the complex changes in hydrology currently observed in the Arctic. Precipitation increases have been observed at many northern sites. However, previous research shows that precipitation increase is not the primary driver for runoff or baseflow increase (Walvoord and Striegl 2007, St Jacques and Sauchyn 2009, Connon et al 2014) and cannot explain contradictory patterns observed in nearby lakes and ponds (e.g. filling/draining lakes in Necsoiu et al 2013), which are typically attributed to subsurface drainage (e.g. Yoshikawa and Hinzman 2003, Avis et al 2011, Jepsen et al 2013). Additional factors such as evapotranspiration and land cover changes may also have non-negligible roles. For example, Connon et al (2014) found that changes in land cover interact with permafrost and winter groundwater discharge, leading to increased mean annual flow in a Canadian subarctic basin. Other processes that are not taken into account in our model, such as cryosuction and geomorphic changes due to the thawing of ice-rich permafrost, may also have an effect on groundwater-surface water interactions.

Changing the permeability, the surface slope and the recharge rates affects the magnitude of the spatial shift in groundwater discharge, but not the timing (figures 6(a)–(c)). In contrast, changing snow density (i.e. winter insulation) alters the timing but not the magnitude of GWDDR (figure 6(d)). The characteristic pattern common to all nine models is due to the formation of the lateral talik and to the development of the perennial supra-permafrost aquifer (figure 2). Thus, while the magnitude and the timing of the shift from upslope to downslope and from AFS to BFS may change based on model parameters, the spatiotemporal patterns are broadly similar across all parameter combinations indicating that the results presented herein are applicable to permafrost settings with a wide range of physical characteristics. Additional results that compare the sensitivity analysis simulations are presented in the supplementary information.

We selected these parameters for sensitivity analysis as our experience and previous work indicates they can significantly impact model outcomes, although other parameters also can affect results. We also tested sensitivity to the size of the model domain, which does not affect the patterns found in this study. Some assumptions, such as a homogeneous, two-dimensional domain, have also been made in this study to simplify the problem and develop a process-based understanding of changes in groundwater discharge in response to climate warming. For example, subsurface heterogeneity, nonlinear land surface slopes, more variable warming, and/or snowfall trends may delay or advance the groundwater discharge shifts simulated in this study.