Introduction

Arctic environments are warming at a faster rate than any other region on earth. This warming is causing concentrated, rapid hydrological changes such as increased freshwater discharge, earlier spring peak flows, increased precipitation and thawing permafrost (Walvoord & Kurylyk, 2016). Climate change influences almost every characteristic of an Arctic watershed snow and rainfall precipitation distributions; vegetation coverage; groundwater storage; permafrost and thawing depth (Hinzman, Bettez, et al., 2005). Permafrost is defined as soil, rock or other natural material that has been frozen for two or more consecutive years (van Everdingen & Association, 1998). We refer to soils that are frozen for less than two years or that thaw every summer, as seasonally frozen soils. Permafrost has been reported to redirect the flow of groundwater (Hinzman, Johnson, Kane, Farris, & Light, 2000). In this context, the role of permafrost connecting and disconnecting groundwater flow and river flow is especially vulnerable to climate change. Apart from permafrost, seasonally frozen soil will likely impact these connections as well. Understanding the effects of seasonally thawing soils on river discharge could provide valuable insights into long-term changes across transitions from permafrost to non-permafrost regions. Moreover as 55-60% of the land surface in the northern hemisphere is currently frozen during winter (Niu & Yang, 2006), with 23.9% of the exposed land area underlain with permafrost (T. Zhang, Barry, Knowles, Heginbottom, & Brown, 2008), understanding how seasonally frozen soils affect the flow of water is crucial to predict hydrological responses to a changed climate.
Based on a synthesis of multiple arctic terrestrial studies examining arctic freshwater processes across different hydrophysiographical regions, Bring et al., (2016) concluded that warming-induced increases in the active layer thickness will likely lead to changes in the storage capacity of groundwater thereby altering river flow dynamics. Walvoord & Kurylyk, (2016) reviewed multiple arctic water flow models and showed that while groundwater exchange and subsurface connectivity is predicted to increase locally, these models are inconclusive on how surface connectivity and river flow dynamics will change at basin scales. Hinzman, Yoshikawa, & Kane, (2005) suggest that evapotranspiration in the Arctic will increase as temperatures increase, leading to dryer soils and lower river flows. Prowse et al., (2015) put forward that the ecological transition, from tundra to boreal, is strongly hydrologically mediated. All of these studies show there is no single dominant permafrost thaw effect on the hydrologic cycle; instead several interacting changes that exacerbate other changes to create complex responses, which differ by region and are typically hard to predict. Predicting such complex and interacting changes in the Arctic is one of the key challenges in hydrology (Peel & Blöschl, 2011; Tetzlaff, Carey, McNamara, Laudon, & Soulsby, 2017). Therefore, observation based approaches are crucial to reveal the ongoing change trajectories, identify dominant processes and build reliable models that can project change trajectories into the future. Although still considered sparse, river discharge is the most commonly available observation throughout the Arctic (Laudon & A. Sponseller, 2017) and contains integrated signals of watershed processes affected by Arctic warming (Bring et al., 2016). In this study, we aim to quantify long-term trends in how watersheds release water (i.e. trends in watershed storage-discharge relationships) throughout Northern Sweden and evaluate if these trends can be attributed to changes in the spatial extent and timing of frozen soils.
Several studies have previously investigated long-term trends in storage-discharge relationships in the Arctic (Bogaart, Van Der Velde, Lyon, & Dekker, 2016; Brutsaert, 2008; Lyon & Destouni, 2010; Lyon et al., 2009; Sjöberg, Frampton, & Lyon, 2013; Watson, Kooi, & Bense, 2013). Typically, these studies assumed linear storage-discharge relationships during winter baseflow conditions and found that groundwater flows more easily (i.e. resistance to flow reduces) into rivers with more pronounced Arctic warming. Lyon et al. (2009) related such changes to an increased active aquifer depth of between 0.7-1.3 cm/y in Northern Sweden. Still, under non-base flow conditions (i.e. following a rainfall or snowmelt event) the relationship between river discharge and water storage is typically nonlinear (Brutsaert & Nieber, 1977; Kirchner, 2009; Wittenberg & Sivapalan, 1999), and cannot easily be related to active flow depths (Bense, Ferguson, & Kooi, 2009). However, under wetter conditions even more pronounced effects of frozen soils on river discharge are expected as seasonal frost hampers infiltration of melt and rainwater into the deeper groundwater, and impacts groundwater outflow into rivers (Ploum et al, 2019; Walvoord & Kurylyk, 2016). Frozen soils are expected to seasonally alter the hydrological connectivity within watersheds by redirecting water dominantly through shallow (above the frozen layer) and deep flow routes (below the frozen layer) towards rivers (Ploum et al., 2019). A change in hydrological connectivity typically alters the functional form of the storage-discharge relationship (i.e. degree of non-linearity) (Lyon & Destouni, 2010).
How does the hydrologic response of Arctic and sub-Arctic watersheds change as the climate warms and the extent of frozen soils recedes? Following up on the study of Ploum et al., (2019), we expect that under thawed conditions deep groundwater, shallow groundwater and overland flow paths all contribute to discharge, while under frozen top soil conditions, shallow flow paths are dominant, although deep groundwater can still contribute. This increase in flow path diversity is expected to occur over both the long-term as permafrost thaws and summer active layer increases with a warming climate, as well as on a seasonal timescale when seasonal frozen soil thaw starts earlier as spring occurs earlier. Previous studies have shown when catchments become wet and the diversity of flow routes increases, increasingly non-linear storage discharge relations are observed (Brutsaert & Nieber, 1977).
Based on these studies, our hypothesis is that as seasonally frozen soils thaw and recede, flow path diversity and thus hydrologic connectivity increases, thereby increasing the non-linearity of the storage-discharge relationship. In this current study, our objective is to test and expand on this hypothesis by quantifying trends and spatio-temporal differences in storage-discharge relationships for sixteen watersheds within Northern Sweden throughout the years of 1950 and 2018. Northern Sweden, has had strong temperatures increases from the start of the 1800’s to the 2000’s of almost 0.1°C/decade, with a cooling period occurring between 1940’s to 1970’s (Klingbjer & Moberg, 2003). Multiple proxies for seasonal and intra-annual differences in extent and depth of frozen soils are used to test whether the observed trends and patterns in storage-discharge relationships can be related to thawing soils.

Methods

Concepts: Recession analysis, storage-discharge relationships, and frozen soils.
Recession analysis is a well-established hydrologic method to examine storage-discharge relationships of watersheds and also offers the advantage of being relatively insensitive to meteorological forcing (Tallaksen, 1995). Recession analysis relates the rate of decline of river discharge to the absolute river discharge (recession curve). Under the assumption of a unique relationship between storage, discharge and a closed water balance, recession curves quantify the storage-discharge relationship (Brutsaert, 2008; Kirchner, 2009). Therefore recession curves can be used to better understand the watershed storage-discharge relationship’s response to climate change, connecting changes in subsurface groundwater flow with changes in surface flows (Ploum et al., 2019; Shaw & Riha, 2012; Wrona et al., 2016).
Recession curves are typically quantified by fitting a nonlinear storage-discharge relationship to hydrograph recessions during periods when changes in discharge reflect changes in catchment water storage (Brutsaert & Nieber, 1977; Kirchner, 2009; Ploum et al., 2019).
\(\frac{\text{dQ}}{\text{dt}}=-Q\frac{\text{dQ}}{\text{dS}}\approx-\alpha Q^{\beta}\)(1)
Here, we extend eq. 1 to alleviate the constrains of no-rain and no-evapotranspiration conditions, which allows us to include and account for periods with small amounts of precipitation and evapotranspiration relative to discharge.
\(\log\left(-\frac{\text{dQ}}{\text{dt}}\right)-\log\left(1+\frac{E}{Q}-\frac{P}{Q}\right)=\log\left(\alpha\right)+\beta\ log(Q)\)(2)
Q is discharge of the river [mm/d], dQ/dt is the rate of discharge decline during the recession [mm/d/d], P is precipitation [mm/d] and E is evapotranspiration [mm/d].dQ/dS is the sensitivity of discharge [d-1] (Kirchner, 2009). In recession analysis, the hydrograph recession is typically plotted in a recession plot with\(\log\left(-\frac{\text{dQ}}{\text{dt}}\right)-\log\left(1+\frac{E}{Q}-\frac{P}{Q}\right)\ \)against\(\ \log\left(Q\right)\) (Equation 2). When Eq.2 is fitted to the data in this log-log plot, \(\alpha\) represents the intercept\({[mm}^{1-\beta}d^{\beta-2}]\), and \(\beta\) is the slope of the recession curve [-]. This technique has been applied to many regions to further understanding of the groundwater and river discharge relationship (Brauer, Teuling, Torfs, & Uijlenhoet, 2013; Dralle, Karst, Charalampous, Veenstra, & Thompson, 2017; Lyon et al., 2015).
Watershed scale recession slopes of rivers can be conceptually interpreted as a measure of hydrologic connectivity as illustrated by a series of buckets (Figure 1). A watershed with one dominant flow path where discharge increases linearly with storage behaves similar to a bucket with a single spigot representing a linear reservoir (\(\beta\) ≈ 1) (Figure 1). Observed examples are watersheds with a deeply incised river during baseflow conditions, a confined aquifer below permafrost, or shallow water flow above a frozen soil (Brutsaert & Hiyama, 2012; Lyon & Destouni, 2010; Ploum et al., 2019). In a drained unconfined aquifer both the pressure as well as the saturated thickness control flow (Troch et al., 2013). Such a system can be represented by a bucket with multiple evenly distributed equally sized spigots (flow paths): flow not only depends on the pressure exerted by storage on the spigots but also by the number of spigots. Such a bucket behaves as a nonlinear reservoir with β = 1.5. A bucket with increasing number of spigots or increasing size of spigots towards the surface will haveβ > 1.5. Examples are found in (Kirchner, 2009), (Karlsen et al., 2019) and (Brauer et al., 2013). A special case is when the resistance of the spigots decreases exponentially towards the surface yielding an exponential reservoir (\(\beta\ \)= 2). This is frequently observed in relatively flat catchments (e.g. Bogaart et al., 2016; Brauer et al., 2013). Reservoirs where the resistance declines hyperbolically towards the surface yield β > 2. Several examples are found in literature (Brutsaert & Nieber, 1977; Kirchner, 2009; Troch et al., 2013).
Figure 1
Recession analysis has been used to examine the rate of thawing permafrost, (Brutsaert, 2008; Lyon & Destouni, 2010; Lyon et al., 2009). Brutsaert (2008) focused on baseflow winter (i.e. frozen) conditions when the recession slope (β) can be assumed 1 (linear reservoir). Lyon & Destouni, 2010 used coefficient α (assuming β = 1) to infer aquifer thickness during late summer, when the active layer extend is maximum in permafrost soils. However, under spring and early summer conditions, we expect that the spatial heterogeneity of actively thickness and hydrological connectivity between aquifers challenges the assumption of linear reservoir behavior (β = 1). Therefore, we focus on temporal changes in recession slope (β ) during spring and summer, which we relate to changes in hydrological connectivity using the interpretation of Figure 1. The spring and summer periods are later separated in order to untangle effects of seasonal soil frost and permafrost on river recessions, with spring period having potential seasonal soil frost and summer as the period with the greatest active layer and without/less seasonal soil frost.
Study Sites and data
Our study sites are situated in Northern Sweden. The sixteen watersheds were chosen because of expected presence of regions with permafrost in the past and present (Brown et al., 1997; Gisnås et al., 2017; Zhang et al., 1999), widespread occurrence of seasonal frozen soils and no current or past known obstruction of the waterways by human intervention following Sjöberg et al. (2013).