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

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:

$$\begin{align} \frac{\partial \Theta}{\partial t} = \frac{\partial}{\partial z} \left[ K(\Theta)\left( \frac{\partial \left(\Psi(\Theta) -z \right)}{\partial z}\right) \right] \end{align} \tag{ 1 }$$

where $\Theta$ is volumetric soil water content (m³ m⁻³), t is time (s), K is hydraulic conductivity (m s⁻¹), Ψ(Θ) 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)

$$\begin{align} \Psi(\Theta) = \frac{1}{\alpha}\left[\theta^{\frac{-1}{m}} - 1~\right]^{\frac{1}{n}} \end{align} \tag{ 2 }$$

where α (m⁻¹) and n (−) are the van Genuchten parameters and $m = 1-\frac{1}{n}$, θ is relative saturation (−) defined as

$$\begin{align} \theta = \frac{\Theta - \Theta_\mathrm r}{\Theta_\mathrm s - \Theta_\mathrm r}. \end{align} \tag{ 3 }$$

$\Theta_\mathrm s$ is the volumetric soil water content at saturation (m³ m⁻³) and $\Theta_\mathrm r$ is the residual water content (m³ m⁻³) at K(Θ) = 0. K(Θ) is parameterized by

$$\begin{align} K(\Theta) = K_\mathrm s \theta^{\frac{1}{2}}\left[1 - \left(1 - \theta^{\frac{1}{m}}\right)^m\right]^2 \end{align} \tag{ 4 }$$

where $K_\mathrm s$ is saturated hydraulic conductivity (m s⁻¹). 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:

$$\begin{align} S_0~\frac{\partial h}{\partial t} = \left(\frac{\partial }{\partial x},\frac{\partial}{\partial y}\right) \cdot \left( k(h)\left(\frac{\partial }{\partial x},\frac{\partial}{\partial y}\right) h \right) \end{align} \tag{ 5 }$$

where S₀ is the specific storage coefficient (m³ m⁻³), h is hydraulic head (ground water table elevation (m) for unconfined conditions), k is transmissivity (m² s⁻¹). This right hand side consists of divergence $\left(\frac{\partial }{\partial x},\frac{\partial}{\partial y}\right) \cdot$ and gradient $\left(\frac{\partial }{\partial x},\frac{\partial}{\partial y}\right)$ 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:

$$\begin{align} C_\mathrm{eff}(\Theta , T_\mathrm{s})\frac{\partial T_\mathrm{s}}{\partial t} = \frac{\partial}{\partial z}\left( \lambda (\Theta , T_\mathrm s) \frac{\partial T_\mathrm s}{\partial z}\right) \end{align} \tag{ 6 }$$

where $C_\mathrm{eff}(\Theta , T_\mathrm{s})$ is effective heat capacity (J m⁻³ K⁻¹) which includes phase change energy (Daanen and Nieber 2009), $T_\mathrm{s}$ is soil temperature (K) and $\lambda(\Theta , T_\mathrm{s})$ is the effective thermal conductivity (W m⁻¹ K⁻¹). $C_\mathrm{eff}(\Theta , T_\mathrm{s})$ is calculated as the arithmetic mean of volumetric heat capacities of dry soil, liquid water and ice. $\lambda(\Theta , T_\mathrm{s})$ is the arithmetic mean of $\lambda(\Theta , T_\mathrm{s})$ of air and geometric mean of $\lambda(\Theta , T_\mathrm{s})$ of dry soil, liquid water and ice. Here, advection is not considered, however, soil moisture movement can still affect the heat transport. The heat equation is solved for the entire soil column (saturated and unsaturated zones). The upper boundary condition is set by the temperature at the ground surface (calculated by the snow heat transfer scheme) and the lower boundary condition is given as a constant geothermal heat flux. Solution of (6) includes estimation of the liquid water fraction $S_\mathrm{l}(\Theta , T_\mathrm{s})$ which is used to couple heat transfer with vertical and lateral water transport. The former coupling is done by substituting $\theta$ in (4) with $S_\mathrm{l}\cdot\theta$. $S_\mathrm{l}(\Theta , T_\mathrm{s})$ is estimated as a function of van Genuchten parameters, and can be considered as a typical parameterization for 'unfrozen water curve'. Coupling with the lateral moisture transport (5) is accomplished by adjusting the transmissivity k. If soil temperatures are not considered, for an unconfined aquifer, k can be estimated by

$$\begin{align} k(h) = K_{xy}(h-H_\mathrm b) \end{align} \tag{ 7 }$$

where K_xy is the saturated hydraulic conductivity in lateral directions (m s⁻¹) and $H_\mathrm 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

$$\begin{align} k(h) = K_{xy}(h-H_\mathrm{b}-h_\mathrm{f}) \end{align} \tag{ 8 }$$

where $h_\mathrm{f}$ is the thickness of the frozen soil within the aquifer (m). $h_\mathrm{f}$ is defined by

$$\begin{align} h_\mathrm{f} = \int_{H_\mathrm{b}}^{h}\Xi\left( T_\mathrm{s}(S_l \leqslant S_\mathrm{tr})-T_\mathrm s\right) dz \end{align} \tag{ 9 }$$

where Ξ() is a Heaviside step function and $S_\mathrm{tr}$ is a threshold value for liquid soil water fraction ($S_\mathrm 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 multi-century 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⁻¹) are distinguished

$$\begin{align} R = R_\mathrm s + R_\mathrm i + R_\mathrm b \end{align} \tag{ 10 }$$

where R is total runoff, $R_\mathrm s$ is surface runoff (infiltration excess from liquid precipitation and snow melt), $R_\mathrm i$ is interflow (slope-induced seepage in the unsaturated zone) generated in the upper soil column (up to 3 m deep) and $R_\mathrm 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_\mathrm{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_\mathrm{pot}$ is split up into potential transpiration ($E_\mathrm{tr,pot}$) and potential evaporation ($E_\mathrm{ev,pot}$) according to vegetation cover fraction (portion of the cell area covered by vegetation). Actual evapotranspiration ($E_\mathrm{act}$) is calculated as the sum of actual transpiration ($E_\mathrm{tr,act}$) and evaporation ($E_\mathrm{ev,act}$). $E_\mathrm{tr,act}$ is withdrawn from the soil within the root zone (0.4 m deep in this study) according to:

$$\begin{align} E_\mathrm{tr , act} = E_\mathrm{tr , pot} \int_{0}^{z_\mathrm{r}} f_\mathrm{r} (\Psi) d_\mathrm{r} (f_\mathrm{r} , z) dz \end{align} \tag{ 11 }$$

where $z_\mathrm{r}$ is depth of the root zone (m), $f_\mathrm{r}$ is a reduction factor (−) due to soil moisture availability (dryness stress) and $d_\mathrm{r} (f_\mathrm{r} , z)$ is effective root density distribution (−). The numerical implementation of $d_\mathrm{r}$ includes calculation of compensation from lower parts of the root zone due to dryness stress in the upper parts. The difference ($E_\mathrm{tr,un}$) of calculated $E_\mathrm{tr,act}$ and $E_\mathrm{tr,act}$ for which no dryness stress is applied in (11) is included in the calculation of actual evaporation:

$$\begin{align} E_\mathrm{ev , act} = \left( E_\mathrm{ev,pot} + E_\mathrm{tr , un}\right) \int_{0}^{z_e} S_l \theta \left( 1 - \frac{z}{z_e} \right) dz \end{align} \tag{ 12 }$$

where $z_\mathrm e$ is maximum evaporation depth (m). This parameterization assumes that evaporation energy dissipates linearly with depth. Together, precipitation, evapotranspiration and runoff form a watershed-wide water balance:

$$\begin{align} \Delta S = P - E_\mathrm{act} -R \end{align} \tag{ 13 }$$

where ΔS is a change in subsurface water storage (mm d⁻¹). In this study we consider annual means of the water balance and its components.

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²) 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.

Zoom In Zoom Out Reset image size

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).

$K_\mathrm s = K_{xy}$(m s−1) a	1D a		2D b
High (10−5)	Warming	Cooling	Warming	Cooling
Medium (10−6)	Warming	Cooling	Warming	Cooling
Low (10−7)	Warming	Cooling	Warming	Cooling

^a Saturated hydraulic conductivity, see equations (4) and (8). ^b See figure 1.

We perform numerical experiments for three values of saturated hydraulic conductivity ($K_\mathrm 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_\mathrm 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² watershed is divided into 16 grid cells (250 × 250 m). In both cases vertical spacing between layers ranges from 0.08 m (upper soil) to 0.5 m (deep soil). At the lateral boundaries we assume zero flux conditions. In this configuration explicit interactions between river channel and ground water table are possible, i.e. baseflow can be generated in addition to surface runoff and interflow. The geometry of our setup is similar to the idealized slopes and watersheds simulated in other studies (McKenzie and Voss 2013, Evans and Ge 2017, Lamontagne-Hallé et al 2018).

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) 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⁻¹.

Zoom In Zoom Out Reset image size

We chose these data for our atmospheric forcing to ensure that under steady 'cold' climate the 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 (−):

$$\begin{align} c = \begin{cases} \displaystyle\frac{\mu(T_{a,c}) - \mu (T_{a,c} \leqslant 0) + \Delta T }{\mu(T_{a,c}) - \mu (T_{a,c}\leqslant0)} & T_{a,c}\leqslant 0 \\ 1 & T_{a,c} \gt 0 \end{cases} \end{align} \tag{ 14 }$$

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:

$$\begin{align} T_a(t,i) = \begin{cases} \displaystyle\frac{i}{n}T_{a,c}(t)+\frac{n-i}{n}T_{a,w}(t) & \textrm{warming} \\[6pt] \displaystyle\frac{i}{n}T_{a,w}(t)+\frac{n-i}{n}T_{a,c}(t) & \textrm{cooling} \end{cases} \end{align} \tag{ 15 }$$

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.075 $^{\circ}\textrm{C}\ \textrm{yr}^{-1}$) is consistent with the slope of linear trend for average air temperature during the cold period in the continuous permafrost zone (Brown et al 2002) and in the Central Interior Alaska climate division (Bieniek et al 2012) 0.06 ± 0.02 $^{\circ}\textrm{C}\ \textrm{yr}^{-1}$ and 0.07 ± 0.01 $^{\circ}\textrm{C}\ \textrm{yr}^{-1}$ (spatial mean ± spatial standard deviation) respectively (calculated from the CRU dataset (Harris et al 2014) over 1970-2018 time period). Under the steady 'warm' conditions the watershed becomes completely thawed.

For 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⁻². 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⁻³.

Under the warming scenario all three $K_\mathrm s$ cases experience an increase in total runoff (R) concurrent with a decrease in evapotranspiration ($E_\mathrm{act}$). These changes are most pronounced with higher $K_\mathrm s$ (table 2). The new 'warm' steady evapotranspiration is reached when the transition ends regardless of the $K_\mathrm s$ (figure 3). Total runoff reaches new equilibrium values shortly after the transition period (10–15 years). In the high $K_\mathrm s$ case, changes in total runoff are solely due to the increase in interflow ($R_\mathrm i$), while the surface runoff ($R_\mathrm s$) remains the same. As $K_\mathrm s$ decreases, the increase in interflow grows. Surface runoff drops in the medium and in particular the low $K_\mathrm s$ case but remains constant in the high $K_\mathrm s$ case. After the transition surface runoff is the same for all three cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. Water balance and runoff components for cooling and warming scenario and high, medium and low saturated hydraulic conductivity ($K_\mathrm s$) under the 1D catchment configuration (see figure 5 for details). All values are based on the annual means of evapotranspiration (E_act ), runoff (R), interflow ($R_\mathrm{i}$) and surface runoff ($R_\mathrm{s}$).

		Cooling			Warming
Variable		High $K_\mathrm s$ a	Medium $K_\mathrm s$ a	Low $K_\mathrm s$ a	High $K_\mathrm s$ a	Medium $K_\mathrm s$ a	Low $K_\mathrm s$ a
$E_\mathrm{act}$	Δ b (mm yr−1)	25	24	12	−25	−24	−12
	Range c (mm yr−1)	25	24	12	25	25	12
	Trend d (mm yr−2)	−0.6	−0.7	−0.5	0.5	0.7	0.4
R	Δ (mm yr−1)	−25	−24	−12	25	24	12
	Range (mm yr−1)	127	110	28	75	81	40
	Trend (mm yr−2)	2.3	2.4	1.0	−1.4	−1.8	−1.0
Ri	Δ (mm yr−1)	−25	−36	−80	25	36	80
	Range (mm yr−1)	127	121	96	75	81	80
	Trend (mm yr−2)	2.3	2.3	3.4	−1.4	−1.8	−2.3
Rs	Δ (mm yr−1)	0	12	68	0	−12	−68
	Range (mm yr−1)	1	12	70	1	12	70
	Trend (mm yr−2)	0.0	−0.3	−2.4	0.0	0.2	1.7

^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.

For all three $K_\mathrm 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_\mathrm 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_\mathrm 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).

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_\mathrm b$); instead of two in the 1D simulations (only surface and interflow), the transition between two 'steady' regimes have different timescales and short-term 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.

Table 3. Water balance and runoff components for cooling and warming scenario and high, medium and low saturated hydraulic conductivity ($K_\mathrm s$) under the 2D catchment configuration (see figure 5 for details). All values are based on the annual means of evapotranspiration ($E_\mathrm{act}$), runoff (R), interflow ($R_\mathrm{i}$) and surface runoff ($R_\mathrm{s}$).

		Cooling			Warming
Variable		High $K_\mathrm s$ a	Medium $K_\mathrm s$ a	Low $K_\mathrm s$ a	High $K_\mathrm s$ a	Medium $K_\mathrm s$ a	Low $K_\mathrm s$ a
$E_\mathrm{act}$	Δ b (mm yr−1)	25	24	12	−25	−24	−12
	Range c (mm yr−1)	25	24	12	25	25	13
	Trend d (mm yr−2)	−0.6	−0.7	−0.5	0.5	0.6	0.4
R	Δ (mm yr−1)	−25	−24	−12	25	24	12
	Range (mm yr−1)	150	90	67	55	44	30
	Trend (mm yr−2)	−3.3	1.4	0.7	1.7	−0.8	−0.9
$R_\mathrm i$	Δ (mm yr−1)	32	35	6	−32	−35	−6
	Range (mm yr−1)	150	94	47	50	75	32
	Trend (mm yr−2)	−3.5	−2.7	−1.1	1.9	2.8	1.2
$R_\mathrm s$	Δ (mm yr−1)	0	10	69	0	−10	−69
	Range (mm yr−1)	1	10	69	1	10	69
	Trend (mm yr−2)	0.0	−0.3	−2.4	0.0	0.2	1.8
$R_\mathrm b$	Δ (mm yr−1)	−57	−69	−87	57	69	87
	Range (mm yr−1)	60	70	113	60	90	88
	Trend (mm yr−2)	0.3	1.9	3.7	−0.4	−3.2	−3.6

^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.

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_\mathrm 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_\mathrm s$ simulation the peak is preceded by a short-term drop in total runoff. For the other two $K_\mathrm{s}$ cases this drop is small and the peak is mostly due to the rapid increase in baseflow, while in the high $K_\mathrm{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_\mathrm{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_\mathrm{s}$ in warm conditions is caused by the fact that when $K_\mathrm 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_\mathrm s$. This means that the baseflow can be exfiltrated into the river channel within a larger portion of the watershed (compared to higher $K_\mathrm 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_\mathrm 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_\mathrm 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_\mathrm 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.

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 near-surface 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.

Zoom In Zoom Out Reset image size

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_\mathrm 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_\mathrm s$ decreases.

Table 4. Difference in mean annual specific soil moisture storage (mm m⁻²) between 'cold' and 'warm' quasi-steady states for different model configurations and depths.

	1D			2D
Depth (m)	High $K_\mathrm s$ a	Medium $K_\mathrm s$ a	Low $K_\mathrm s$ a	High $K_\mathrm s$ a	Medium $K_\mathrm s$ a	Low $K_\mathrm s$ a
0–0.4	24	30	44	24	30	44
0.4–3	466	435	162	373	179	100
3–70	−965	−925	−411	−4040	−1880	−70
Total	−475	−460	−205	−3643	−1671	74

^a See table 1 for values.

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_\mathrm 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_\mathrm 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_\mathrm 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_\mathrm s$ baseflow is consistent throughout the year.

Zoom In Zoom Out Reset image size

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_\mathrm s$. In the '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_\mathrm 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%.

In terms of the thermal regime under the steady climates we observe no significant difference between different $K_\mathrm 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_\mathrm s$ (even though our model considers only heat conduction and no advection). With increasing $K_\mathrm 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.

Zoom In Zoom Out Reset image size