Changing groundwater discharge dynamics in permafrost regions

Permafrost thaw due to climate warming modifies hydrological processes by increasing hydrological connectivity between aquifers and surface water bodies and increasing groundwater storage. While previous studies have documented arctic river baseflow increases and changing wetland and lake distributions, the hydrogeological processes leading to these changes remain poorly understood. This study uses a coupled heat and groundwater flow numerical model with dynamic freezing and thawing processes and an improved set of boundary conditions to simulate the impacts of climate warming on permafrost distribution and groundwater discharge to surface water bodies. We show a spatial shift in groundwater discharge from upslope to downslope and a temporal shift with increasing groundwater discharge during the winter season due to the formation of a lateral supra-permafrost talik underlying the active layer. These insights into changing patterns of groundwater discharge help explain observed changes in arctic baseflow and wetland patterns and are important for northern water resources and ecosystem management.


Introduction
Arctic and subarctic regions are experiencing major hydrological changes as a result of climate warming (Walvoord and Kurylyk 2016). Over recent decades, pronounced atmospheric warming rates at high latitudes (McBean et al 2005, Turner et al 2007 have contributed to a significant decrease in the distribution and thickness of permafrost (i.e. ground with temperatures perennially below 0 • C) in many northern regions (Walsh et al 2005, Romanovsky et al 2010, and this degradation is expected to accelerate during the next century (Slater andLawrence 2013, Chadburn et al 2017). Permafrost acts as an impermeable barrier which confines groundwater flow to the seasonally-thawed active layer above permafrost, in sub-permafrost aquifers, and through perennially unfrozen zones called taliks (Woo 2012). Thus, as permafrost thaws, groundwater storage and connectivity increases and may impact wetlands and surface water bodies that rely on groundwater exchange (Walvoord and Kurylyk 2016). However, due to the diversity and remoteness of northern environments, little is known about the impacts of permafrost thaw on groundwatersurface water interactions.
Many studies have reported increases in baseflow and mean annual discharge in northern rivers (e.g. Walvoord and Striegl 2007, Déry et al 2009, St Jacques and Sauchyn 2009, Sjöberg et al 2013, Connon et al 2014. Although various explanations for these trends have been proposed, several studies have shown that increases mainly occur when stream discharge is primarily groundwater-derived, suggesting that the changes are related to permafrost thaw (Walvoord and Striegl 2007, St Jacques and Sauchyn 2009, Duan et al 2017. Increased baseflow due to groundwater discharge may be caused by (1) increased hydrogeologic connectivity due to new groundwater pathways formed by permafrost thaw (Walvoord and Striegl 2007, St Jacques and Sauchyn 2009), (2) a local increase of precipitation and atmospheric moisture that enhances groundwater recharge (Connon et al 2014, Duan et al 2017, and/or (3) a temporary decrease in water storage capacity due to changes in land-surface cover (Connon et al 2014).
Increased hydrogeologic connectivity can occur both vertically and laterally. Increased vertical connectivity is caused when supra-and sub-permafrost aquifers connect, such as when an open talik forms below surface water bodies (Bense et al 2012, Walvoord et al 2012. Increased lateral connectivity is caused by the formation of lateral taliks, which are known to exist both within (French 2013) and above (Sloan and van Everdingen 1988, Delisle 2007, Zhang et al 2008 permafrost and function as thin, perennial aquifers (e.g. Sloan andvan Everdingen 1988, Hiyama et al 2013). In general, the influence of permafrost thaw on vertical and lateral hydrogeologic connectivity remains poorly understood.
In contrast to many arctic rivers experiencing increased baseflows, lentic surface water bodies exhibit complex drying and wetting patterns. Net area decrease is observed or predicted for many arctic ponds and wetlands (Yoshikawa and Hinzman 2003, Avis et al 2011, Necsoiu et al 2013 and lakes (Smith et al 2005, Carroll et al 2011, while other regions are undergoing large increases in lake surface area (Korosi et al 2017). Smith et al (2005) observed a decline in Siberian lake cover in discontinuous permafrost regions, but an increase in continuous regions. Conversely, Carroll et al (2011) reported that the highest net loss of lake area in the Canadian Arctic is located in continuous permafrost. Jepsen et al (2013) found that groundwater and active zone dynamics impact lake growth or shrinkage. These differences suggest that the processes underpinning arctic hydrological trends are not well understood and cannot be only attributed to local climate conditions and permafrost distribution.
Several studies have simulated changes in cold regions groundwater discharge in a warming climate through coupled groundwater flow and heat transfer models, and have concluded that groundwater discharge to surface water bodies is expected to increase in most cases (Ge et al 2011, Bense et al 2012, Kurylyk et al 2014a, Evans and Ge 2017, Evans et al 2018. However, as Ge et al (2011) and acknowledge, most previous studies have employed fully saturated conditions and simplified boundary conditions (BC) that can lead to unrealistic groundwater recharge and discharge values. Choosing a combination of BCs that is representative of reality is a persistent challenge in cold regions groundwater modelling (Sanford 2002, Walvoord and. To our knowledge, no long-term modelling studies of permafrost thaw and coupled groundwater dynamics have directly incorporated snow cover insulation effects and seasonal variations of groundwater recharge into a groundwater flow and energy transport model, nor have they investigated the spatial and temporal evolution of groundwater discharge due to permafrost thaw. To fill this knowledge gap, we ask: what are the major processes leading to groundwater discharge changes in permafrost environments and how will those processes influence future baseflow trends? We investigate this question using a heuristic modelling approach, in which theoretical hillslopes in archetypal permafrost environments are used to investigate relationships between permafrost thaw, supra-permafrost aquifer thickness and groundwater discharge. Specifically, we (1) employ improved surface-layer BCs, (2) predict the short-and long-term effects of climate change on the spatial distribution of land surface groundwater discharge due to permafrost thaw, and (3) analyze the evolution of seasonal patterns of baseflow.

Modelling approach and model domain
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, 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 −13 m 2 and exponentially reduces to 10 −19 m 2 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, twodimensional 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 crosssection 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.

Boundary conditions
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 −2 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 −1 and 24.3 mm month −1 , 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.
For the land surface thermal BC, previous modelling studies typically used either a specified ground surface temperature (GST) boundary (Bense et al 2012 or a specified AT with an underlying thermal boundary layer (McKenzie et al 2007, Evans andGe 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.

Modelling sequence
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 −1 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 −1 , 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 −1 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.

Sensitivity analysis
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.

Spatial evolution of groundwater discharge and hydrogeological connectivity
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)).
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.

Changing baseflow seasonality
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 Table 1. Values used for the eight additional simulations. ℎ , , rat , and snow stands for the horizontal and vertical permeabilities, the ratio of precipitation that is set as recharge in the model, and the snow density respectively.

Model parameters
Higher value Default value Lower value Permeability ( ℎ, × or / 10) ℎ = 1 × 10 −12 m 2 = 1 × 10 −13 m 2 ℎ = 1 × 10 −13 m 2 = 1 × 10 −14 m 2 ℎ = 1 × 10 −14 m 2 = 1 × 10 −15 m 2 Surface slope (Slope × or / 2) Slope = 4% Slope = 2% Slope = 1% Recharge rate (R rat ± 50%)   is 9%, which is within the range of arctic precipitation projections (Kattsov et al 2005). 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 suprapermafrost 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 surfacesubsurface 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.
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.

Sensitivity analysis
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, twodimensional 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.

Conclusions
This study advances numerical modeling investigations of groundwater flow in permafrost regions by assessing both the spatial and temporal patterns of groundwater discharge changes and provides an understanding of the mechanisms causing these variations. For this, an improved set of boundary conditions that considers many processes generally ignored in cryohydrogeological models has been applied to an archetype groundwater model that includes heat transport and phase change of pore water.
Permafrost thaw leads to increased hydrogeologic connectivity and increases in river baseflow. These changes are most pronounced when a lateral talik forms and provides a perennial conduit through which groundwater can flow to the downgradient river. Further, increased supra-permafrost thickness (i.e. a lowering of the permafrost table) causes groundwater discharge to migrate downgradient towards basin outlets, such as rivers, rather than into upgradient areas (e.g. wetlands). Although other processes (e.g. evapotranspiration, land cover change) may impact groundwater discharge patterns, these results provide a phenomenological explanation for the observed paradox of arctic landscape wetting (increased baseflow) and drying (shrinking wetlands and lakes).
Recents reports have highlighted the rapidly changing hydrology across the pan-arctic basin (AMAP 2017). Since arctic regions are expected to warm faster than anywhere else in the world, it is important to understand the impacts of climate warming on arctic and subarctic hydrology. Knowing where and when groundwater discharge variations occur is key for predicting the future of arctic river ecosystems, carbon fluxes, and for water management of northern communities that use surface water bodies as their primary water consumption source (Lemieux et al 2016, Walvoord and. For example, Kuhn et al (2018) show the importance of thermokarst ponds as carbon source to the atmosphere. The draining of these ponds may strongly alter the carbon sink capacity of northern wetlands. Terrestrial and aquatic ecosystems in arctic environments also depend on the spatial and temporal distribution of surface and/or subsurface water. This study reveals that the upgradient regions of watersheds will experience surficial drying, which may impact wetland ecosystems. Increased winter groundwater discharge may also increase winter stream temperatures and shift the latitudinal limits of the migratory patterns for certain fish species. Our results predict that the most pronounced changes in arctic and subarctic hydrology and cryohydrogeology are expected to be realized over the next few decades, and thus there is a great need for studies that couple modelling approaches with field investigations in permafrost regions in the near future (e.g. Ghias et al 2017).