Water balance response of permafrost-affected watersheds to changes in air temperatures

Observations show increases in river discharge to the Arctic Ocean especially in winter over the last decades but the physical mechanisms driving these changes are not yet fully understood. We hypothesize that even in the absence of a precipitation increase, permafrost degradation alone can lead to increased annual river runoff. To test this hypothesis we perform 12 millennium-long simulations over an idealized hypothetical watershed (1 km2) using a distributed, physically based water balance model (Water flow and Balance Simulation Model, WaSiM). The model is forced by both a hypothetical warming defined by an air temperature increase of 7.5 ∘C over 100 years, and a corresponding cooling scenario. To assess model sensitivity we vary soil saturated hydraulic conductivity and lateral subsurface flow configuration. Under the warming scenario, changes in subsurface water transport due to ground temperature changes result in a 7%–14% increase in annual runoff accompanied by a 6%–20% decrease in evapotranspiration. The increase in runoff is most pronounced in winter. Hence, the simulations demonstrate that changes in permafrost characteristics due to climate warming and associated changes in evapotranspiration provide a plausible mechanism for the observed runoff increases in Arctic watersheds. In addition, our experiments show that when lateral subsurface moisture transport is not included, as commonly done in global-scale Earth System Models, the equilibrium water balance in response to the warming or cooling is similar to the water balance in simulations where lateral subsurface transport is included. However, the transient changes in water balance components prior to reaching equilibrium differ greatly between the two. For example, for high saturated hydraulic conductivity only when lateral subsurface transport is considered, a period of decreased runoff occurs immediately after the warming. This period is characterized by a positive change in soil moisture storage caused by the soil moisture deficit developed during prior cooling.


Introduction
At regional and pan-Arctic scales observations have shown a significant increase in total annual river discharge to the Arctic Ocean as well as spatial and seasonal changes in river runoff over the last decades (Shiklomanov et al 2013, 2020, Tan and Gan 2015).Observations of combined river discharge from the six largest Russian north flowing rivers have shown an increase of 7% over the period 1936-1999(Peterson et al 2002)), with >70% of the contribution due to higher river flow in winter, although constructed reservoirs in these river basins significantly distort seasonal discharge and complicate our understanding of river flow due to climate variability (Shiklomanov andLammers 2009, Stuefer et al 2011).With effects of reservoir regulation removed, even greater annual runoff increases of 9.4 km 3 yr −1 (13% over  were found (Shiklomanov et al 2013).Recent estimates show that the increase in river flow from Eurasia to the Arctic Ocean has continued into the 21st century with the new maximum recorded flow in 2007 (Holmes et al 2016).River flow has also increased in North America at a rate of 0.9 km 3 yr −2 (7% over  for Mackenzie, Yukon, Peel and Back (Shiklomanov and Lammers 2011) and by 8.4 km 3 yr −1 (18% over  for all rivers in northern Canada (Déry et al 2016).The most recent assessment of observed river flow to the Arctic Ocean from Eurasia and North America shows a 5.1 km 3 yr −1 or 9% increase in total influx over 1964-2015(Déry et al 2016, Holmes et al 2018).
Atmospheric moisture transport and associated precipitation patterns play a major role in runoff generation, and its changes can significantly alter the hydrological regime of rivers flowing to the Arctic Ocean.However, aggregated over the Arctic and large river basins, annual precipitation, which is typically the most important water balance component for runoff generation, does not show a substantial change to support the observed increasing trend in annual river flow (Berezovskaya et al 2004, Adam and Lettenmaier 2008, Bring and Destouni 2011, Bring et al 2016).The most consistent and significant increase in river flow throughout the Eurasian pan-Arctic is observed during the cold season (Shiklomanov et al 2013, Tananaev et al 2016).Similar tendencies in winter discharge were found in northern North America (Déry et al 2016).Some increase in snowfall during fall and early winter in Siberia (Wegmann et al 2015) and in the Canadian Arctic (Kopec et al 2016) associated with sea ice decline can partly explain increases in spring river flows but not an increase in winter flow.In contrast to precipitation, the increase in air temperature across the Arctic has been widely and consistently documented (Biskaborn et al 2019).The air temperature rise leads to significant changes in the regional cryosphere with corresponding impacts on the hydrological regime.For example, the watershed of the Little Chena Creek, a tributary to the Yukon River through the Tanana River, with a long runoff record has experienced a significant increase in runoff over the last 40 years.The increase is two times higher than the increase in annual precipitation for that period, but the watershed has also experienced a significant change in air temperatures which is most pronounced in the winter (figure S1 in supplementary materials (available online at stacks.iop.org/ERL/16/084054/mmedia)).
Permafrost, defined as a layer of soil that stays below 0 • C for more than two consecutive years (Muller 1945, Permafrost Subcommittee 1988) underlays approximately 24% of the exposed land in the Northern Hemisphere (Zhang et al 1999, Brown et al 2002).The recent rise in air temperature in the Arctic and sub-Arctic regions resulted in widespread thawing and degradation of permafrost (Romanovsky et al 2017).Many local and regional studies have highlighted the effect these changes in permafrost can have on catchment-scale water balance (Woo et al 2008, Connon et al 2014, Tananaev et al 2016, Walvoord and Kurylyk 2016) providing a plausible explanation for increases in runoff as winter air temperatures rise.Permafrost acts as an impermeable hydrological barrier reducing soil water storage capacity and constraining subsurface flow (Woo et al 2008, Walvoord and Kurylyk 2016, Tananaev et al 2020) thus exerting a major control on runoff processes and the interaction of surface water and sub-surface water storage (Rawlins et al 2019).However, a clear quantitative assessment of the role of permafrost thaw on the observed increases in river flow is still missing.Comprehensive reviews of recent advances in permafrost hydrology are given by Woo et al (2008), Walvoord and Kurylyk (2016), Fabre et al (2017) and Tananaev et al (2020).To translate the conceptual understanding of how coupled subsurface heat and moisture transfer influences the terrestrial hydrological cycle into quantitative estimates of this influence, many modeling studies have been performed in recent years (e.g.Frampton et al 2013, McKenzie and Voss 2013, Evans and Ge 2017, Jafarov et al 2018, Lamontagne-Hallé et al 2018).However, most of these studies do not include interactive nearsurface hydrological processes (evapotranspiration, surface runoff, interflow generation) or they resolve only moisture transport within the saturated zone.
Model intercomparison studies suggest that on a pan-Arctic scale models exhibit considerable difficulties in realistically representing hydrological processes (Gädeke et al 2020).In addition, projections of the future response of Arctic watersheds to climate change show that though on average the models predict Arctic drying, they diverge considerably in their projected trajectories of soil moisture, evapotranspiration and runoff (Andresen et al 2020).The large spread in model projections can be attributed to a myriad of factors including differences in parameterizations of hydrological processes, model structure, spatial representation of soil parameters, and forcing data.A detailed examination of isolated individual processes on a small watershed scale provides a way forward to improve representation of permafrostaffected watersheds on a global scale.
In this study, we perform controlled numerical experiments with a detailed, physically-based hydrological model (Schulla 1997) in a hypothetical watershed to investigate in a controlled way how surface water balance is altered by a transition from frozen to thawed soil conditions and vice versa.We hypothesize that even in the absence of increased precipitation, air temperature-driven permafrost degradation can lead to significant river flow increases due to enhanced subsurface moisture transport.This, in turn, leads to a drier near-surface soil causing a decrease in evapotranspiration.We hypothesize that this process could be acting over large regional scales and could have a significant contribution to the overall river flow increases observed in recent years in the Arctic.To explore this hypothesis and understand the sensitivity of the results we perform model experiments with varying soil saturated hydraulic conductivity.In addition, we investigate the differences in simulations with and without the temperature-dependent lateral subsurface moisture transport within the saturated zone to understand the impact of neglecting lateral temperature-dependent subsurface water transport in modern Earth System Models.
The controlled model experiments are designed to isolate the effects of permafrost degradation and aggradation on the watershed's water balance.We aim to achieve this by assuming changes in air temperature solely in the cold (winter) period within the year.This allows us to keep potential evapotranspiration and liquid/solid precipitation partitioning constant within the entire simulation period.We designed our experimental watershed setup to be as generic as possible.We look at an idealized basin in an idealized continental environment with typical precipitation that represents higher elevation watersheds in northern North America and Eurasia.We acknowledge that we cannot cover all the variety of environment types and variables impacting runoff, and the suggested experimental design can be furthermore expanded to assess the long-term sensitivity of permafrost-affected watersheds to the variety of parameters under changing climate.Nevertheless, our goal for the experimental design in this study is to provide a basic framework for such assessment.

Model description
We use the Water flow and balance Simulation Model (WaSiM), a well-established and widely used tool for simulating the spatial and temporal variability of hydrological processes in complex river basins WaSiM (Schulla 1997).The model is distributed, deterministic, and physically based.The governing equations and their numerical approximations are described fully in Schulla (2019).Here, we describe those processes most relevant for the scope of this study.For each cell heat transfer within the entire soil column and moisture transfer within the unsaturated zone are modeled in one dimension vertically.Lateral surface and subsurface (in the saturated zone) water transport between the cells are modeled separately.The vertical movement of water in the unsaturated part of the soil column is described by the Richards equation: where Θ is volumetric soil water content (m 3 m −3 ), t is time (s), K is hydraulic conductivity (m s −1 ), Ψ(Θ) is soil water potential (m), z is depth below the ground surface (m).The function Ψ(Θ) is modeled using the van Genuchten-Mualem parameterization (Van Genuchten 1980) where α (m −1 ) and n (−) are the van Genuchten parameters and m = 1 − 1 n , θ is relative saturation (−) defined as Θ s is the volumetric soil water content at saturation (m 3 m −3 ) and Θ r is the residual water content (m 3 m −3 ) at K(Θ) = 0. K(Θ) is parameterized by where K s is saturated hydraulic conductivity (m s −1 ).
The upper boundary condition is provided by the vertical moisture flux from the surface (from liquid precipitation or snow melt).The lower boundary condition comes from iterative balancing of vertical fluxes between the solution of the Richards equation and the solution for the lateral flow in the saturated zone which is modeled by a numerical solution of Darcy's law for an uppermost unconfined aquifer: where S 0 is the specific storage coefficient (m 3 m −3 ), h is hydraulic head (ground water table elevation (m) for unconfined conditions), k is transmissivity (m 2 s −1 ).This right hand side consists of divergence • and gradient in lateral Cartesian coordinates (x, y).Boundary conditions (lateral) are assumed to be zero flux in this study, but are adjustable in the model.
The heat transfer in soil which is crucial for permafrost regions is implemented by numerically solving the one-dimensional heat equation in the vertical direction: where where K xy is the saturated hydraulic conductivity in lateral directions (m s −1 ) and H b the elevation of the bottom of the soil column.Generally, frozen ground largely inhibits any water movement since only a small fraction of the water content remains liquid.Therefore, we modify (7) to where h f is the thickness of the frozen soil within the aquifer (m).h f is defined by where Ξ() is a Heaviside step function and S tr is a threshold value for liquid soil water fraction (S l , usually ≈0.9).This modification provides a reduction of the transmissivity when the aquifer is partially frozen.
If the entire aquifer is frozen little to no lateral water movement will occur.
WaSiM provides an array of parameterizations of the processes involved in surface water balance that provides top boundary conditions for the Richards and heat equations and deal with surface energy and water balance.In this study, since we consider multicentury time-scales, we employ simplified approaches for runoff generation and routing as well as for evapotranspiration.Snow heat transport is modeled with the n-factor approach (Kade et al 2006), and snow melt is computed using a simple temperature index approach.For runoff generation three runoff components (mm d −1 ) are distinguished where R is total runoff, R s is surface runoff (infiltration excess from liquid precipitation and snow melt), R i is interflow (slope-induced seepage in the unsaturated zone) generated in the upper soil column (up to 3 m deep) and R b is baseflow (exfiltration of ground water into the rivers).The depth of the interflow withdrawal is roughly determined by the field capacity (matric potential of ≈3.45 m).By definition interflow is generated from soil to the surface.Hence, any excess water from layers below ≈3.45 m can not be withdrawn as interflow since it will not be able to reach the surface (given the matric potential requirement).
Baseflow is generated in any cell if the groundwater table is above the subgrid channel bed.It is assumed that the subgrid channel has an open talik and baseflow can be generated even if the soil in the actual cell is frozen.Subgid channel and saturated zone interactions are controlled by the colmation factor which in this study is assumed to be equal to the saturated hydraulic conductivity of the soil.All three components of runoff from each cell are then routed separately by a kinematic wave approach (inverse Manning formula) through subgrid river channels.
Evapotranspiration is derived from potential evapotranspiration (E pot ) based on the Hamon parameterization (Hamon 1961).This approach is chosen since it only depends on air temperature and does not require additional meteorological variables.E pot is split up into potential transpiration (E tr,pot ) and potential evaporation (E ev,pot ) according to vegetation cover fraction (portion of the cell area covered by vegetation).Actual evapotranspiration (E act ) is calculated as the sum of actual transpiration (E tr,act ) and evaporation (E ev,act ).E tr,act is withdrawn from the soil within the root zone (0.4 m deep in this study) according to: where z r is depth of the root zone (m), f r is a reduction factor (−) due to soil moisture availability (dryness stress) and d r (f r , z) is effective root density distribution (−).The numerical implementation of d r includes calculation of compensation from lower parts of the root zone due to dryness stress in the upper parts.The difference (E tr,un ) of calculated E tr,act and E tr,act for which no dryness stress is applied in ( 11) is included in the calculation of actual evaporation: where z e is maximum evaporation depth (m).This parameterization assumes that evaporation energy  4) and ( 8).
b See figure 1.
dissipates linearly with depth.Together, precipitation, evapotranspiration and runoff form a watershedwide water balance: where ∆S is a change in subsurface water storage (mm d −1 ).In this study we consider annual means of the water balance and its components.

Model experiments
To investigate how permafrost degradation and aggradation fundamentally influences the water balance of small permafrost-affected watersheds we perform 12 idealized numerical experiments with WaSiM varying soil saturated hydraulic conductivity, lateral subsurface flow and atmospheric forcing (table 1) .The experiments are applied to a hypothetical watershed defined by a uniform slope of 1 • over a length of 4000 m and a uniform width of 250 m (1 km 2 ) with a single subgrid river channel along the slope (figure 1).All experiments assume homogeneous soil texture with no anisotropy for the entire 70 m deep soil column and the same vegetation type (tussock tundra) and associated model parameters.
We perform numerical experiments for three values of saturated hydraulic conductivity (K s ).These are chosen to encompass the range of drainage conditions typically encountered in continuous permafrost soils in the Arctic.For all experiments we consider isotropic conditions (K s = K xy ).For each of these cases we run two experiments: (a) we consider a 1D case with only vertical subsurface water transfer, i.e. the entire watershed is represented in the model by one single grid cell (figure 1).Hence, no interactions between the river channel and groundwater table are possible within the watershed; therefore runoff is only generated as surface runoff and interflow.These experiments allow us to estimate the watershed response in the same manner as most of Land and Earth System Models do, which include heat and moisture transfer in soil only in the vertical direction.(b) We consider a 2D case where subsurface water transfer is possible in both the vertical and horizontal direction along the slope (figure 1).The 1 km 2 watershed is divided into 16 grid cells (250 × 250 m).In both cases vertical spacing between layers ranges The model is forced by near-surface air temperature (T a ) and precipitation (P) data.For each of the six model configurations described above we apply two atmospheric forcing scenarios: a warming scenario leading to permafrost degradation and a cooling scenario causing permafrost aggradation.Forcing scenarios are constructed as linear transitions (figure 2) of air temperature over 100 years between two steady climates, referred to as a 'cold' climate and a 'warm' climate (figure 2).The 'cold' climate represents current conditions and consists of annual cycles of air temperature and precipitation that are derived from averaging spatial means of the Central Interior Alaska climate division (over period 1970-2018)  We chose these data for our atmospheric forcing to ensure that under steady 'cold' climate the (current) and 'warm' climate.Only winter temperatures are altered between the 'cold' and 'warm' climate while the annual precipitation cycle is assumed identical for both climates.Time series of air temperature (bottom) transitioning over a 100 year period from a steady 'cold' climate to a steady 'warm' climate.Also shown are winter (defined as season when air temperatures are negative), summer (positive air temperatures) and annual means.The case of the transition from steady 'warm' climate to a steady 'cold' climate is not shown but derived analogously.hypothetical watershed is completely underlain by permafrost (mean annual ground surface temperature is −1.7 • C) while the watershed is not too cold so that the transition to the completely thawed state occurs within a reasonable time for the entire soil column of the watershed.The resulting monthly annual cycles (figure 2) are interpolated to the model time step (30 min) with the first four Fourier coefficients to produce an annual cycle and guarantee its smoothness and lack of noise.Then, a linear adjustment is applied to ensure a smooth transition from one year to the next i.e. the beginning and the end of two consecutive annual cycles have the same T a and P.
Then, we construct a 'warm' climate annual cycle where we ensure that the increase in air temperature does not affect the summer period (defined as the period with positive daily air temperatures).Onset and duration of summer as well as the summer air temperature are the same for any given year in the simulations(figure 2) and only winter temperatures are being changed.Thus the snow melt starts at the same time, annual, snow and rain precipitation amount is preserved, and potential evaptranspiration are the same for both the 'warm' and the 'cold' climates.The annual cycle of the 'warm' climate under this assumption is then obtained by multiplying the half-hourly air temperatures of the 'cold' climate with the following coefficient c (−): where µ(T a,c ) denotes the 'cold' climate's mean annual air temperature, and ∆T the mean temperature difference between 'cold' and 'warm' climates.The annual cycle of precipitation are equal for both the 'cold' and 'warm' climate and our treatment of air temperatures preserves the individual liquid/solid parts.We apply ∆ T = 7.5 • C to ensure that permafrost is absent in most of the soil column (mean annual ground surface temperature is 2.8 • C) in the 'warm' climate.This allows us to achieve a complete transition between the thawed and frozen state of the entire domain while no other aspects of the annual water and heat balances are affected.
In principle, there are alternative ways of forcing the model so that similar changes in the surface temperature and resulting transition from frozen to thawed state are achieved, e.g. by changing snow depth.However, in contrast to our approach, such changes in forcing would lead to additional changes in the water balance.These unrelated changes would no longer allow us to isolate the effect of atmospheric warming on the water balance (through temperature-dependent soil moisture transport mechanism).
We construct warming and cooling scenarios defined by a linear transition between 'cold' and 'warm' steady climates and vise versa (figure 2).The transition spans n = 100 years and air temperature (T a (t, i)) at any given time t : t ∈ [0, 1] (normalized time within a calendar year) and year i : i = 1, 2, 3, …, n within the transition is derived as: where T a,w (t) and T a,c (t) are near surface air temperatures in the 'warm' and 'cold' climates, respectively.The length of the transition period is chosen to be 100 years in order to ensure a reasonable time for calculations and a sufficiently small year to year change in air temperature during the transition.The choice of ∆T (equation ( 14)) and a resulting rate of change during the transition period (0.075For all six model configurations (regardless of the forcing scenario) the model is initialized at the start of the hydrologic year (1 October) with a uniform soil temperature distribution equal to the mean annual ground surface temperature for the 'warm' climate (2.8 • C) and a ground water table at 10 m depth (soil column above is completely dry).The bottom boundary condition for heat transfer is set to a constant heat flux of 0.085 W m −2 .Lateral and bottom boundary conditions for soil moisture transport are set to zero water flux.The model is then spun up for 1000 years under 'warm' climate conditions to reach an equilibrium state.After the spin up the model is forced with a cooling scenario.The end state of the cooling scenario is then used as the initial condition for the warming scenario.Both scenarios consist of 100 years of steady climate, 100 years of transition, and 900-1400 years of the opposite steady climate.The length of the simulation depends on when the relative difference in the decadal averages of soil temperature and moisture distributions of two consecutive decades is no greater than 10 −3 .

Water balance components 1D cases
Under the warming scenario all three K s cases experience an increase in total runoff (R) concurrent with a decrease in evapotranspiration (E act ).These changes are most pronounced with higher K s (table 2).The new 'warm' steady evapotranspiration is reached when the transition ends regardless of the K s (figure 3).Total runoff reaches new equilibrium values shortly after the transition period (10-15 years).In the high K s case, changes in total runoff are solely due to the increase in interflow (R i ), while the surface runoff (R s ) remains the same.As K s decreases, the increase in interflow grows.Surface runoff drops in the medium and in particular the low K s case but remains constant in the high K s case.After the transition surface runoff is the same for all three cases.
For all three K s cases, total runoff shows a temporary spike roughly 95 years in the transition.This peak is caused by a rapid increase in interflow coinciding with a rapid increase in liquid soil moisture (through talik development) at 0.4-3 m depth (figure 4 and figures S2-S4 in supplementary materials).Soil moisture in the upper 0.4 m of the soil (where surface runoff and evapotranspiration are withdrawn) decreases smoothly throughout the whole transition period.While the moisture content for the upper most layer (up to 0.08 m) does not change in the case of high K s , for the other two cases it decreases significantly within the transition period facilitating a decrease in surface runoff (figure S2 in supplementary materials).
Under the cooling scenario (transition from 'warm' climate to 'cold' climate, see figure 2) all variables return to the initial values of the warming scenario (table 2), however the steady-state values of evaptranspiration and total runoff are generally reached faster, and in the cases with lower K s well before the end of the transition period.In addition, a rapid drop in total runoff and a subsequent increase can be observed in the first 30 years after the onset of the transition (figure 3).The underlying reason for this behavior is the rapid drop in interflow.This is a direct consequence of the fact that freezing starts from the top.As soon as an initial perennially frozen layer is established the underlying soil within the interflow generation layer (3 m deep) freezes rapidly while the freezing front propagation slows down with depth.In the deeper layers where interflow is generated, the soil moisture perennially freezes reducing the volume of total liquid moisture available for the interflow generation to only the seasonally thawed layer.On the other hand, the near-surface soil layers where the moisture is spent on evaptranspiration and surface runoff the soil becomes wetter throughout the year (figure 4).    1) in 1D configuration.Liquid and frozen components are stacked.Shaded area shows transition period between steady climates with vertical dashed line marking the year when mean annual air temperature crosses 0 • C.

Water balance components in 2D case
Overall, the experiments with increased model complexity (2D simulation, figure 5) show similar behavior to the corresponding 1D simulations (figure 3).The 'warm' and 'cold' steady water balances have similar differences for evapotranspiration and total runoff (tables 2 and 3).However, since in the 2D simulations the runoff generation is represented with three mechanisms: surface, interflow and baseflow (R b ); instead of two in the 1D simulations (only surface and interflow), the transition between two 'steady' regimes have different timescales and shortterm magnitudes.The evolution of spatially averaged (along the slope) soil moisture, temperature and liquid fraction are similar to those of the 1D cases and can be found in figures S5-S7 in supplementary materials.
In the warming scenarios evapotranspiration in the 2D cases follows the same patterns as in the 1D cases in terms of the long term change and the transient changes and are closer in magnitude for all K s cases.In contrast, runoff exhibits a significantly different transient behavior (figure 5).The differences in timing of the short-term peak in runoff are greater in the 2D case.In the high K s simulation the peak is preceded by a short-term drop in total runoff.For the other two K s cases this drop is small and the peak is mostly due to the rapid increase in baseflow, while in the high K s case the baseflow gradually increases over 1200 years after the transition to the 'warm' climate.In addition, the overall length of the period when the surface water balance is not at equilibrium (non-zero change in subsurface storage) is longer (up to three times) than for the experiments with the 1D configuration.As the heatwave from increased air temperature propagates into the deeper soil it allows for soil moisture to engage in interflow generation more effectively similar to the 1D configuration.However, as soon as the saturated zone begins to thaw the lateral subsurface transport becomes more effective due to the k(T) relationship ( 8).The interflow generation decreases and the baseflow starts to increase in total runoff.In the most extreme case (high K s ) this takes over 900 years (while the change in storage stays at zero).This switch between interflow and baseflow occurs due to the increased downward soil moisture flux in the upper part of the slope since the water in the saturated zone is now allowed to move laterally following the hydraulic gradient in the unsaturated zone.This increased downward flux leaves the upper 3 m of soil where interflow can be generated with less soil moisture.In the lower part of the watershed the extra water provided by lateral movement is exfiltrated into the river channel thus generating baseflow.Surface runoff behaves similarly to the experiments with 1D configuration.Increased baseflow under lower K s in warm conditions is caused by the fact that when K s is low-the vadose zone is shallow (due to the relatively restricted water movement).This restriction comes from the lower colmation associated with low K s .This means that the baseflow can be exfiltrated into the river channel within a larger portion of the watershed (compared to higher K s cases) and interflow generation is limited to only the upper part of the watershed where the vadose zone is deeper than 3 m.
Under the cooling scenario we observe a similar behavior where long-term changes in evapotranspiration and surface runoff in the 2D configuration are similar to those in the 1D configuration.In contrast, total runoff, though showing similar overall long-term change between steady climates, exhibits a significantly different transient behavior than in 1D configuration, especially for the case of high K s .Instead of a sharp drop in total runoff, we observe a slight drop followed by a significant increase that later gradually recesses to the steady 'cold' value.
For the two other K s cases the total runoff behavior differs from 1D configuration experiments in terms of the magnitude of the initial drop in total runoff (tables 2 and 3).The transient behavior in the high K s case is due to the long smooth recession of the baseflow and a sharp increase in the interflow near the time when mean annual air temperature crosses 0 • C.These changes are caused by the establishment of the newly formed frozen layer that reduces the influx from the soil surface to the deeper soil.However, since the deeper layers are gradually freezing from the top, soil moisture in the previously thawed unconfined aquifer underneath the freezing front is still allowed to move laterally.The aquifer slowly freezes while its unfrozen part is still involved in lateral subsurface transport until the latter stops and baseflow can no longer be generated.

Equilibrium water balance regimes under steady conditions
Total annual evapotranspiration and total runoff for steady 'warm' and 'cold' climates in the 2D experiments differ by less than 0.5% from the respective values in 1D experiments.In both, 1D and 2D simulations the water balance regime in the 'warm' climate is characterized by higher annual runoff and lower evapotranspiration compared to the 'cold' climate due to better connectivity between near-surface soil layers and deeper soil.The differences range from ≈6% to ≈20% for the evapotranspiration and ≈7% to ≈14% for the total runoff (table 2 and  figures 3 and 5).The difference in the connectivity is also apparent in the annual cycle of the nearsurface soil moisture (figure 6).Under the 'cold' regime, due to overall higher moisture content in the upper 0.4 m of soil, the total volume of liquid water available for evapotranspiration is higher.In addition the soil between 0.4 and 3 m is mostly frozen so that interflow generation at these levels is suppressed.Under the 'warm' regime the total soil moisture is lower in the upper 40 cm and the amount of available water for evapotranspiration during summer is reduced.Although the average annual water content for the depth between 0.4 and 3 m is generally lower under the 'warm' regime, this moisture is perennially liquid and is always available for interflow generation.
With regard to changes to total soil moisture storage significant differences between 1D and 2D cases can be observed (table 4).Though the upper soil layers (up to 3 m deep) have similar differences between 'cold' and 'warm' regimes for both 1D and 2D cases, the deeper soil layers (below 3 m) have significantly different values between 1D and 2D cases.This distinction between 1D and 2D configuration is most prominent in the case of high K s .This is due to significantly lower water storage in the 'cold' state, since after the establishment of the impermeable frozen layer during the cooling scenario the moisture stored in the lower soil column is still allowed to move and generate runoff (through lateral subsurface transport as baseflow) before the soil column is mostly frozen.This effect decreases in magnitude as K s decreases.
The difference in the subsurface connectivity between the 'cold' and 'warm' regimes is also reflected in the annual runoff cycles (figure 7).Almost the entirety of the surface runoff contribution to the annual water balance for all model configurations, K s cases and forcing scenarios occurs in the beginning of the warm period and consists mostly of snow melt water.Under the 'cold' climate only medium and high K s have significant summer and early fall runoff that consists entirely of interflow and no runoff is present in winter month for all cases.Under the 'warm' regime on the other hand, the runoff persists during the winter for all cases and the snow melt peak of surface runoff is significantly lower.The difference between 1D and 2D model configurations is most apparent in the medium K s case where towards the end of summer interflow is replaced by baseflow and regains its dominant position until the end of the winter.In the case of high and low K s baseflow is consistent throughout the year.
In terms of runoff components, the 'warm' regime for 1D cases is dominated by interflow where it constitutes ≈60% to ≈70% of the total runoff with it is contribution increasing with K s .In the In terms of the thermal regime under the steady climates we observe no significant difference between different K s cases and model configurations (figure 8 and figures S3 and S6 in supplementary materials).However, due to the differences in in annual cycle of the soil moisture (figure 6) the effectiveness of the heat transfer in the upper soil column changes with the K s (even though our model considers only heat conduction and no advection).With increasing K s and decreasing soil moisture the ratio between summer and winter thermal conductivity also increase both under the 'cold' and 'warm' steady climates.In addition, in all cases, under the 'cold' climate this ratio is lower which results in a larger thermal offset (difference between the ground surface temperature and ground temperature at the bottom of seasonal freezing/thawing).This increased thermal offset allows for permafrost temperature at the depth of zero annual amplitude to be lower by 1.3 • C than ground surface temperature and also contributes to less time needed for the watershed to reach steady state under the cooling scenario.

Discussion
In our numerical experiments we isolate the mechanism of temperature dependent moisture transport which, we hypothesize, can alone be a sufficient mechanism that allows for higher runoff under warmer climate conditions.The isolation of this mechanism is achieved by keeping summer air temperatures and annual precipitation constant while only the winter temperatures are allowed to change.Other studies that investigated the long-term reaction of subsurface flow to the increase in air temperature (e.  temperature throughout the year.For the timescales that are considered in this and other studies only the total yearly increase in additional heat at the surface is relevant.However, other studies utilized constant or periodic functions of pressure heads and/or recharge flux at the soil-atmosphere interface (consistent with our assumption of constant annual precipitation).This approach does not consider moisture-dependent evapotranspiration.Thus, in contrast to those studies, our experimental setup allows us to evaluate the watershed's response to air temperature increase in greater detail: not only in terms of groundwater discharge but considering all water balance components (13) and runoff partitioning (10).
The similarity of the long-term water balance response of the simulated watersheds in 1D and 2D model configurations suggests that the application of simpler and less computationally intensive models that do not consider lateral subsurface transport can produce comparable estimates of the changes in water balance due to changes in climate over millennial time scales assuming that the overall change in air temperature is sufficient to transition a watershed between mostly frozen and mostly thawed states.However, the differences in the watershed's behavior during the transition period is drastically different for 1D and 2D cases.This suggests that 1D coupled soil heat and moisture transfer schemes often used in Earth System Models and Land Models (Andresen et al 2020) might be sufficient to represent differences in water balance between equilibrium states but will most likely incorrectly predict transient changes in water balance.
In cases with lateral subsurface transport (2D configuration) after the initial perennially frozen layer is established, the liquid water in the new subpermafrost aquifer is still allowed to move.If the 'warm' equilibrium water table elevation is low enough, this subpermafrost aquifer becomes essentially an unconfined aquifer with zero additional soil moisture flux from the top boundary (figure 9) after the onset of freezing.If the local river channel along the watershed is large enough to maintain a through talik underneath it-the subpermafrost aquifer is free to drain through this talik into the river channel.Such conditions are consistent with our setup where a subgrid river channel is assumed to have interaction with the underlying soil regardless of its temperature state.Since the elevation of the water table at the top of the watershed cannot be maintained by the incoming infiltration it begins to drop while the subpermafrost aquifer is drained into the river channel.Depending on the K s value a substantial deficit of soil moisture storage can be developed in the higher elevation part of the watershed before it freezes up completely.During the subsequent warming, this deficit needs to be refilled before the 'warm' equilibrium water table can be restored.Inclusion of a parameterization for cryosuction (Unold and Derk 2017) in the model could reduce this deficit in soil moisture storage, however, the effect of cryosuction diminishes with the increase in K s (Watanabe et al 2011) and lower soil moisture in the unsaturated zone (Watanabe et al 2012) following a faster lowering of the groundwater table.
In contrast to other modeling studies (e.g.Frampton et al 2013, Evans and Ge 2017, Lamontagne-Hallé et al 2018), our experiments suggest a different baseflow pattern under warming scenarios since, unlike in the other studies, a considerable portion of the watershed can be unsaturated due to the differences in model initialization and watershed geometry (figure 9).In this study, we do not consider Yedoma (late Pleistocene ice-and organic-rich silty sediments) or other massive ground ice formations that are typically found on the flat terraces and floodplains of the major Arctic rivers (Strauss et al 2013) where excessive subsurface water storage is apparent.Instead, our experiments can be interpreted as typical conditions for the regions of the permafrost-affected river basins at high elevations (above regional water tables).In these regions unsaturated soil (consistent with storage deficit observed in this study) can be found in epigenetic (Permafrost Subcommittee 1988) permafrost cores (e.g.French and Shur 2010, Stephani et al 2014, Cable et al 2018).However, the issue of the ground ice distribution across the Arctic and associated water content remains highly uncertain (Hugelius et al 2014), and does not allow us to assess how typical are permafrost conditions exhibiting soil moisture deficit presented in this study.
In our experiments, we consider three values of K s while all the other soil parameters are assumed identical for all cases based on the values typical for sandy loam.However, in real soils, a plethora of highly covariant hydraulic parameters vary between different soil types, not just saturated hydraulic conductivity.This assumption contributes to an increased interflow with the increase of K s under the 'warm' steady state.Since other parameters (e.g.drainage density, recession constant, van Genuchten parameters) are kept constant between the cases they allow higher sensitivity for the interflow (with respect to K s ) than for baseflow outside of their characteristic case (medium K s ).In addition, lower K s values mean less intensive exchange between the river network and the saturated zone leading to an increase in water table elevation upstream.This increase results in a larger area for interaction between the saturated zone and river network which increases baseflow and decreases interflow (shallower vadose zone at the hilltop).However, since in this study the primary interest is in the near-surface process changes (surface runoff and evapotranspiration) and the subsurface (both interflow and baseflow), this specific effect in runoff fractioning (increased interflow in lower K s cases) does not diminish the validity of the main conclusions about the total water balance and soil moisture storage.Nevertheless, these inconsistencies warrant further investigation on the overall sensitivity of such watersheds towards individual parameters and their combinations.
Similar to the other long-term modeling studies (e.g.Voss 2013).In addition, in WaSiM it is not restricted to only one unconfined aquifer, however, it is not possible to dynamically separate into several aquifers within the modeled domain (models implementing full 3D schemes for moisture and heat transport avoid this limitation).This might lead to no explicit lateral flow in the suprapermafrost aquifer in the cooling scenario in our simulations, although this is somewhat compensated by the interflow generation above the newly established permafrost table.
In terms of the millennial scale increase in baseflow our results agree well with those from similar modeling studies (e.g.Evans and Ge 2017, Lamontagne-Hallé et al 2018), regardless of the observed difference in the subsurface water storage under 'cold' climate with this study.In terms of the total runoff and water balance components our results also agree with the conceptual understanding (e.g.Walvoord andKurylyk 2016, Tananaev et al 2020).In addition, it must be noted that the general picture (steady regimes) is qualitatively consistent with observations [figure 7

Conclusions
In this study we demonstrate how, through coupled heat and moisture subsurface transfer, Arctic watersheds react to changes in air temperature.The increase in air temperature and subsequent permafrost degradation results in the increased moisture transport from the surface to the deeper soil.The results of our simulations suggest that this process redistributes soil moisture to the deeper soil leaving the near surface layers dryer than under the warmer climate.This results in increased annual runoff (especially in the winter) and decreased evapotranspiration.When the temperature trend is reversed the connectivity between near surface and deeper soil greatly reduces.This reduction leads to an increase in soil moisture in the near surface layers and an increase in evapotranspiration.The deeper soil deprived of the moisture flux from the layers above releases the excess moisture while it gradually freezes until a new equilibrium is reached.Thus our simulations prove that over long timescales increase in runoff and decrease in evapotranspiration can be caused solely by the increase subsurface moisture transport during warming and the reverse changes in water balance occur under cooling.
The results of our simulations suggest that a simplified approach (common in Land and Earth System Models) that does not include temperaturedependent lateral subsurface flow is able to produce differences in water balance components between equilibrium states similar to those simulated with lateral subsurface flow.However, the transient behavior of the watershed is not adequately captured by the simplified approach.The differences between the approaches become more apparent with the increase in K s of the simulated watershed.The latter has been shown to be one of the most sensitive parameters for the groundwater discharge response to the warming climate (Lamontagne-Hallé et al 2018).
Our experiments with the watershed being initialized from the'warm' equilibrium state and subsequent cooling allows us to demonstrate that typically the 'cold' equilibrium regime is characterized by the deficit in soil moisture storage in the deeper soil layers which is proportional to the soil K s .This indicates that watersheds near the hilltops (above the regional water table elevation) of the Arctic river basins might experience a period of decreased runoff and positive change in water storage during warming.However, in the long-term, the runoff under warming in these watersheds will eventually exceed the initial 'cold' runoff values.
In our simulations we only consider a range of K s values to illustrate the natural variability of the water balance changes.However, a range of other factors and their combinations may significantly influence the timescales and magnitude of the changes in water balance, although the general pattern will most likely remain the same.These factors include the magnitude and rate of the air temperature change in the forcing scenarios for the simulations, topographic attributes of the watershed (slope, profile curvature), heterogeneity of the soil saturated hydraulic conductivity and its anisotropy as well as values of other soil hydraulic parameters and their distribution within a watershed.Further research should assess not only the individual sensitivity to a particular factor but also investigate their joint influence on the changes in water balance driven by the air temperature change.

Figure 1 .
Figure 1.Two model configurations for the hypothetical watershed.1D corresponds to a single cell configuration with no lateral subsurface flow.2D corresponds to 16 cell configuration with subsurface lateral flow.
from the Climate Research Unit (CRU) dataset (Harris et al 2014).Mean annual air temperature for the steady 'cold' climate is −4.8 • C and annual precipitation is 340 mm yr −1 .

Figure 2 .
Figure 2. Model atmospheric forcing.Annual cycles (top) of near surface air temperature and precipitation during the 'cold'(current) and 'warm' climate.Only winter temperatures are altered between the 'cold' and 'warm' climate while the annual precipitation cycle is assumed identical for both climates.Time series of air temperature (bottom) transitioning over a 100 year period from a steady 'cold' climate to a steady 'warm' climate.Also shown are winter (defined as season when air temperatures are negative), summer (positive air temperatures) and annual means.The case of the transition from steady 'warm' climate to a steady 'cold' climate is not shown but derived analogously.

Figure 3 .
Figure 3. Water balance components for 1D simulations.Ta is air temperature (shown are mean annual and mean winter temperatures defined as all days with subfreezing temperatures), Eact is annual actual evapotranspiration, R is total annual runoff, R i is annual interflow and Rs is annual surface runoff.Shaded area marks the transition period from a warm to a cold climate (left panels) and vice versa (right panels) and dashed vertical line marks the year when mean annual air temperature crosses 0 • C. Results are shown for three saturated hydraulic conductivity (Ks) values (table 1).

Figure 4 .
Figure 4. Mean annual volumetric soil moisture content (Θ) averaged for 0 to 0.4 m (evapotranspiration, surface runoff) and 0.4 m to 3 m (interflow, baseflow) for warming and cooling scenarios and three saturated hydraulic conductivity (Ks) values (table1) in 1D configuration.Liquid and frozen components are stacked.Shaded area shows transition period between steady climates with vertical dashed line marking the year when mean annual air temperature crosses 0 • C.

Figure 5 .
Figure5.Water Balance components for 2D simulations.Ta is air temperature (shown are mean annual and mean winter temperatures defined as all days with subfreezing temperatures), Eact is annual actual evapotranspiration, R is total annual runoff, R i is annual baseflow, R i is annual interflow and Rs is annual surface runoff.Shaded area marks the transition period from a warm to a cold climate (left panels) and vice versa (right panels) and dashed vertical line marks the year when mean annual air temperature crosses 0 • C. Results are shown for three saturated hydraulic conductivity (Ks) values (table1).
g. Frampton et al 2013, McKenzie and Voss 2013, Evans and Ge 2017, Lamontagne-Hallé et al 2018) have utilized a constant increase in air

Figure 9 .
Figure 9. Conceptualized sketch of the evolution of the distribution of frozen/unfrozen soil and unsaturated/saturated zone in the idealized watershed in response to the transition between the warm and cold steady climate and vice versa.
and figure S1 in supplementary materials, as well as, e.g.Connon et al 2014, Niu et al 2016, Tananaev et al 2016, Evans et al 2020, King et al 2020] although real watersheds are affected by more contributing factors (precipitation, summer temperatures change, anisotropic and heterogeneous soils, etc) than considered here.

Table 1 .
Numerical experiments performed in the idealized watershed.'Warming' and 'Cooling' refer to an atmospheric warming and cooling scenario, respectively (see text for details).
a Saturated hydraulic conductivity, see equations (

Table 2 .
Water balance and runoff components for cooling and warming scenario and high, medium and low saturated hydraulic conductivity (Ks) under the 1D catchment configuration (see figure5for details).All values are based on the annual means of evapotranspiration (Eact), runoff (R), interflow (R i ) and surface runoff (Rs).
a See table 1 for values.b Difference between the end and the start of the simulation.c Difference between the maximum and the minimum value during the simulation.d Maximum absolute value of the 30 year running trend over the simulation period.

Table 3 .
Water balance and runoff components for cooling and warming scenario and high, medium and low saturated hydraulic conductivity (Ks) under the 2D catchment configuration (see figure5for details).All values are based on the annual means of evapotranspiration (Eact), runoff (R), interflow (R i ) and surface runoff (Rs).
a See table 1 for values.b Difference between the end and the start of the simulation.c Difference between the maximum and the minimum value during the simulation.d Maximum absolute value of the 30 year running trend over the simulation period.

Table 4 .
Difference in mean annual specific soil moisture storage (mm m −2 ) between 'cold' and 'warm' quasi-steady states for different model configurations and depths.Annual runoff cycles at steady climates.Runoff components are stacked.Results for three saturated hydraulic conductivity (Ks, table 1) values and the 1D and 2D model configuration (figure1).Mean annual ground temperature profiles with annual minima and maxima (top), ratio of the summer (λs) and winter (λw) thermal conductivity (middle) and mean annual soil moisture (Θ, includes liquid water and ice) profiles with annual minima and maxima (bottom) for 'cold' and 'warm' steady climates.Results for three saturated hydraulic conductivity (Ks) values (table1) in the 1D model configuration.
'cold' climate interflow contribution in both 1D and 2D (no baseflow occurs in the 'cold' climate) configurations is decreased.Most significantly in the low K s case where interflow contributes to only ≈10% of the total runoff.Under the 'warm' climate in 2D configuration interflow represents ≈10% to ≈40% of the total runoff and baseflow represents ≈30% to ≈60%. Figure 8.
Frampton et al 2013, McKenzie and Voss 2013, Evans and Ge 2017, Lamontagne-Hallé et al 2018) we do not consider dynamic vegetation cover in our numerical experiments.Changes in vegetation due to air temperature and soil moisture can significantly alter the transpiration demand (Raz-Yaseef et al 2017, Zhang et al 2018, Sabater et al 2020) as well as the surface heat and moisture fluxes (Liston et al 2002, Thompson et al 2004, Jafarov et al 2018).However, generally, vegetation succession is mostly dependent on the summer air temperature (Epstein et al 2013, Pearson et al 2013) while in this study we only consider changes in winter temperatures.Unlike other studies that employ models (e.g.McKenzie and Voss 2013, Evans and Ge 2017, Lamontagne-Hallé et al 2018) to simulate the behavior of small permafrost-affected watersheds on millennial timescales the model employed in this study includes representation of near surface processes (e.g.evapotranspiration, dynamic snow melt).However, other models (e.g.Jafarov et al 2018, Lamontagne-Hallé et al 2018) include more detailed representations of the subsurface processes.For instance, the absence of advective heat transfer in WaSiM may lead to longer transition times between 'frozen' and 'thawed' states of the watershed in warming cases (McKenzie and