Modeling the role of preferential snow accumulation in through talik development and hillslope groundwater flow in a transitional permafrost landscape

Elchin E Jafarov; Ethan T Coon; Dylan R Harp; Cathy J Wilson; Scott L Painter; Adam L Atchley; Vladimir E Romanovsky

doi:10.1088/1748-9326/aadd30

1. Introduction

Shifting climate conditions in the Arctic are driving evolution of permafrost landscapes that include changes in, and complex interactions between, plant productivity and species distribution, permafrost degradation, hydrologic processes, and the terrestrial carbon cycle and energy balance (AMAP 2017). Increasing migration of shrubs into tundra could lead to positive changes in surface energy balance and net loss of carbon from the ecosystem (Bret-Harte et al 2002, Cable et al 2016). Shrub-snow interaction could boost winter heating by up to 70% for exposed shrubs, leading to warmer soil temperature below shrubs over winter periods (Sturm et al 2005). A better understanding of the Arctic tundra system under current climate conditions, specifically the role of shrub-snow interaction on talik formation, will ultimately lead to improved projections of future changes in the terrestrial Arctic. In this work, we study the effect of talik on continuous permafrost layer along the Arctic hillslope and its ability to transition continuous permafrost to discontinuous permafrost.

Taliks are perennially thawed ground within permafrost (Jorgenson et al 2010); they stimulate soil carbon respiration and introduce changes to subsurface hydrology (Harden et al 2012, McKenzie and Voss 2013). Knowing how the movement of water changes during talik formation will improve understanding of evolving subsurface thermal regimes, solute transport, and hydrological connectivity (Frampton and Destouni 2015, Sjöberg et al 2016, Ala-Aho et al 2018). Through taliks extend completely through the permafrost and effectively connect the sub-permafrost zone to surface and active layer hydrology. Here we explore fundamental interactions between snowpack depth, subsurface thermal regime, and current climate conditions that drive localized through talik formation on hilly landscapes. Along hilly landscapes, changes in permafrost from continuous to discontinuous may drive significant changes in hydrologic pathways of groundwater movement and the spatial redistribution of water, carbon and nutrients (Walvoord and Kurylyk 2016).

There are several mechanisms besides a warming climate that contribute to talik formation, including heat from standing water bodies (Rowland et al 2011), land disturbance such as wildfires (Jafarov et al 2013) or infrastructure development (Nelson et al 2001), and non-uniform snow depth due to preferential snow distribution (Pomeroy et al 2006, Burn et al 2009). In this study, we explore the effect of non-uniform snow depth along a hillslope on talik development. Snow insulates the ground from colder temperatures during winter and influences the thermal dynamics of the active layer and underlying permafrost (Shiklomanov and Nelson 1999, Zhang 2005, Jafarov et al 2014). Localized areas that consistently accumulate relatively deeper snowpack than their surrounding environment typically have a thicker active layer or may develop a talik (Burn et al 2009, Atchley et al 2016). For example, the use of snow fences to trap windblown snow to protect nearby infrastructure usually results in significant warming of the soil beneath the accumulated snow (Hinkel and Hurd 2006, O'Neill and Burn 2017). Alternatively, areas with tall alder and willow shrub patches, often distributed along streams and across hillslopes, trap snow and are a dominant influence on the spatial distribution of snow depth across the landscape, resulting in the development of taliks (Pomeroy et al 2006).

The boundary between zones of continuous and discontinuous permafrost (Jorgenson et al 2008) is an area of the interplay between vegetation and snow can lead to a non-homogeneous spatial distribution of permafrost. One such area is the Seward Peninsula of Alaska, which has a mixture of tundra type vegetation, shrubs, and trees, classifying this area as a Low Arctic Zone (Bliss 1997). Non-uniform snow depth due to shrub/tundra distribution contributes to the complex spatial heterogeneity in permafrost presence along the landscape at the Seward Peninsula, where talik formation under shrub patches might be one of the dominant mechanisms responsible for creating permafrost discontinuities. To understand how a preferentially distributed snowpack can lead to the formation of taliks, we ask the following questions: can taliks form under the current-day climate and permafrost conditions along Arctic hillslopes, similar to those on the Seward Peninsula? Under what permafrost conditions will through taliks form along hillslopes? What are the consequences of through talik formation on hillslope thermal hydrology?

In this study, we use a recently developed model of permafrost dynamics to answer the questions posed above. There are many models that simulate talik development as a result of deeper snow depth in 1D (McGuire et al 2016, Nicolsky et al 2017). However, 1D models cannot fully describe the hydrothermal changes within a hillslope system. To study how a through talik might change hillslope hydrology, activate connections between shallow and deep flow pathways and affect groundwater storage, we use a 2D hillslope model driven by preferential snow distribution and current-day climatology.

The advanced terrestrial simulator (ATS, v0.86)⁴ is a fully 3D-capable, coupled hydrologic flow and heat transport simulator that includes the soil physics needed to capture permafrost dynamics, including ice, gas, and liquid saturation, multi-layered soil physics, and flow of unfrozen water in the presence of phase change (Coon et al 2016, Painter et al 2016). The ATS model allows heterogeneous soil layering and has been well tested against measured data through a growing history of applications to Arctic permafrost. In particular, Harp et al (2016) showed the importance of subsurface characterizations on projections of permafrost conditions in a warming climate; Atchley et al (2015, 2016) tested sensitivities of hydrothermal parameters in the model, and analyzed the effect of environmental conditions (peat thickness, snow depth, and degree of water saturation) on the active layer thickness in ATS; Sjöberg et al (2016) used ATS to study the role of advective transport of energy in fens; Jan et al (2017) used a mixed-dimensional model structure to model surface/subsurface thermal hydrology in low-relief permafrost regions at a watershed scale; and Schuh et al (2017) utilized ATS to investigate the effect of soil moisture retention curves on silt loam sediment texture and their impact on permafrost hydrothermal conditions.

In this work, we present results of modeling talik formation along a hillslope transect using ATS v0.86. Physically-based boundary conditions are applied to isolate the influence of snow distribution under multiple subsurface thermal regimes that result in talik formation under current climate conditions. First, we briefly introduce the model and the meteorological data used in this modeling exercise, then we outline necessary steps required for the hillslope domain setup. We describe numerical experiments capturing the formation of taliks under a variety of initial permafrost conditions and describe the resulting impacts on the hydrologic and energetic state of the permafrost-rich hillslope. This allows us to evaluate how changes in water movement occur as hydrologic connectivity is established during talik formation.

2. Methods

2.1. ATS model

The ATS is a collection of representations of physical processes (thermal, hydrological, transport, deformation, and so on) designed to work within a flexibly configured modeling framework, Amanzi-ATS. This modeling framework combines Amanzi's infrastructure for solving coupled physical systems on unstructured meshes (Moulton et al 2012) with a multiphysics process model management system called Arcos (Coon et al 2016).

In this study, we configure ATS as described in detail by Painter et al (2016). This configuration of ATS is a distributed, integrated hydrologic model including a surface energy balance, enabling the simulation of snowpack evolution, surface runoff, and subsurface flow in the presence of phase change. The entire system is mass and energy conservative, and is driven by meteorological data and appropriate initial and boundary conditions.

In the surface energy balance (Atchley et al 2015), a balance of incoming and outgoing radiation, latent and sensible heat exchange, and diffusion of energy to the soil are used to determine a snowpack temperature, and therefore predict the evolution of the snowpack. This evolution includes snow aging, density changes, and depth hoar, and delivers fluxes of water and energy to the surface.

Surficial overland flow is solved through the use of the diffusion wave equation based on Manning's surficial roughness coefficient, with modifications to account for frozen surface water. Additionally, transport of energy on the surface is tracked through an advection-dominated transport equation that allows phase change. This surface system is used as a nonlinear boundary condition along the top of the subsurface domain, resulting in continuity of (liquid) water pressure, temperature, and flux of both mass and energy between above- and below-ground processes.

Subsurface thermal hydrology is represented in ATS by a modified Richards' equation coupled to an energy transport equation, which solves for temperature and variably-saturated, water dynamics with phase change. Unfrozen water at temperatures below 0 °C, an important physical phenomenon for accurate modeling of permafrost temperatures (e.g. Romanovsky and Osterkamp 2000), is represented through a constitutive relation (Painter and Karra 2014) based on the Clapeyron equation, which partitions water into the three phases. Density variation with phase change, cryosuction of liquid water into the freezing front in unsaturated soils, and reduction of flow due to freezing water (i.e. the 'freezing equals drying' approximation) are included. Karra et al (2014) demonstrated that this formulation is able to reproduce several variably-saturated, freezing soil laboratory experiments without resorting to empirical soil freezing curves or impedance factors in the relative permeability model.

Other key constitutive equations for this model include a water retention curve, here the van Genuchten–Mualem formulation (Mualem 1976, van Genuchten 1980), and a model for thermal conductivity of the air–water–ice–soil mixture. The method for calculating thermal conductivities of the air–water–ice–soil mixture is based on interpolating between saturated frozen, saturated unfrozen, and fully dry states (Painter et al 2016) where the thermal conductivities of each end-member states is determined by the thermal conductivity of the components (soil grains and air, water, or liquid) weighted by the relative abundance of each component in the cell (Johansen 1977, Peters-Lidard et al 1998, Atchley et al 2015, appendix A). Standard empirical fits are used for the internal energy of each component of the air–water–ice–soil mixture. In particular, the permafrost state is sensitive to choices of thermal conductivity, and parameters needed to fully specify the model are described in detail in table 1; notation is identical to Painter et al (2016).

Table 1. Thermal and hydraulic properties of three subsurface layers used in the model.

Physical properties	Organic matter	Mineral	Cobble-rich mineral soil
Porosity (m³ m⁻³)	0.55	0.45	0.25
Permeability (m²⁾	1 × 10⁻¹²	2 × 10⁻¹³	1 × 10⁻¹⁵
Thermal conductivity, (saturated unfrozen) (W m⁻¹ K⁻¹⁾	0.67	1	1
Thermal conductivity, (dry, unfrozen) (W m⁻¹ K⁻¹⁾	0.07	0.29	0.33
Volumetric heat capacity (10⁶) (J m⁻¹ K⁻¹)	0.4	1.5	1.5
van Genuchten, α (10⁻⁴)	2.93	5.45	5.45
van Genuchten, m	0.212	0.191	0.191

Flexible configuration of the ATS code allows the setup of different numerical discretization schemes and linear and nonlinear solvers. We use a two-point flux discretization for both the surface and subsurface equations; this choice, while less accurate than some, is strictly monotonicity-preserving, which is critical for accurate and efficient solution of freezing soil dynamics. A Nonlinear Krylov Acceleration method (Calef et al 2013), modified to include significant globalization including backtracking and physics-based strategies, was used to solve the resulting globally-implicit coupled system of equations.

2.2. Site and meteorological data descriptions

This work assists interpretation of field research that is being undertaken as part of the Next Generation Ecosystem Experiment, NGEE-Arctic project⁵ . The project aims to improve the representation of Arctic ecosystem processes in Earth System Models by identifying and characterizing important system interactions, poorly understood processes, and landscape heterogeneity. In particular, NGEE-Arctic researchers are performing field investigations across a gradient of continuous and discontinuous permafrost watersheds on the Seward Peninsula, Alaska. The idealized hillslope modeled in this study serves as a generalization of the snow-vegetation interaction observed at the Teller watershed (figure 1), which through geophysical and in situ observations, has been shown to have strong variation in permafrost depth and structure that appears to be correlated with the distribution of tall shrub patches (Dafflon et al 2016).

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** (A) The transect along the hillslope at the Teller site marked by red line. (B) The location of the Teller study site at the Seward Peninsula in Alaska marked by yellow rectangle.
Download figure:
Standard image High-resolution image

Observations suggest the presence of permafrost close to the ground surface at the upslope and downslope portions of the watershed and variable distribution of a deeper permafrost table (greater than 1.5 m), or possibly no permafrost, along much of the midslope portion of the watershed. The organic layer thickness (OLT) varies from the top to the bottom of the watershed, with thick peat and moss layers at the top and the bottom and little to no moss or peat within the midslope under tall willow shrub patches. The low thermal conductivity of the peat layer makes it an effective insulator, decreasing the heat exchange between the atmosphere and underlying soil layers (Jafarov and Schaefer 2016). The collapsed peat plateau wetland area at the bottom of the watershed is characterized by grass and sedge vegetation and thick peat. Snow depth surveys conducted along the watershed indicate that areas with tall willows (1–2 m in height) trap more snow than low stature vegetation, i.e. tundra, grass, or sedge. This conceptualization is shown in figure 2.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Schematics of the modeling transect.
Download figure:
Standard image High-resolution image

The Teller site inspired the conceptual model domain which aims to explore general interactions between vegetation-induced snow depth, ground temperatures, and subsurface hydrology along the typical Arctic hillslope. In this study, we forced the model with detrended meteorological data derived from the Daymet v3 gridded to a 1 by 1 km product (Thornton et al 2017) from 1980–2015 for the region near Teller, AK (figure 1). We compared air temperatures and precipitation from Daymet dataset for Teller with the same data from National Ocean and Atmospheric Administration's (NOAA's) meteorological station at Nome. The comparison shown that Daymet has on average 2.4 °C colder mean annual temperatures and more precipitation. We chose Daymet data because it comes with shortwave radiation, humidity, and wind which are not available through NOAA's station and because we are consider our modeling domain to be broadly representative of hillslopes in transitional permafrost, not a specific site. The processed data includes air temperature (K), incoming shortwave radiation (W m⁻²), rainfall (m s⁻¹), snowfall in snow water equivalent (SWE) (m SWE s⁻¹), relative humidity (-), and wind speed (m s⁻¹). The incoming long-wave radiation is calculated internally as a function of air temperature and relative humidity (Atchley et al 2015).

Figure 3(A) shows averaged yearly statistics for processed air temperatures from the Daymet v3 data⁶ . The mean annual air temperature (MAAT) averaged over the 36 year record is −5.7 °C, and has a positive trend of 0.06 °C per year. Figures 3(B) and (C) show yearly cumulative rain and snow precipitation near Teller, Alaska. The yearly cumulative rain precipitation also shows a positive trend with a rate of 8.4 mm per year. The total average rain precipitation over the 36 year record is approximately 444 mm. The yearly cumulative snow precipitation showed a negative trend of −4.8 mm SWE per year, where the total average snow precipitation over the 36 year record is 333 mm. The linear trends in the data were removed to generate a statistically consistent, typical 'current-day' climatology. It is important to note that the Daymet data do not exactly correspond to local meteorological data collected at the Teller site. We use the Daymet data instead of the local data due to short time series and incompleteness of measured meteorological data at the site.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** The red line with dots corresponds to (A) the yearly averaged air temperature, (B) yearly cumulative rain, and (C) snow precipitation from 1980–2015 from the Daymet v3 dataset corresponding to the Teller site, Alaska. The blue, solid lines correspond to best-fit linear trends, and the black dashed lines correspond to the mean values.
Download figure:
Standard image High-resolution image

2.3. Mesh and material properties

We designed the mesh and chose material properties to represent a generalized low Arctic hillslope partially based on observations from the Teller watershed. The conceptual model is a hillslope transect with tundra vegetation along the top and bottom, and shrubs along the middle. The upslope end of the hillslope is considered a groundwater divide along the top of a ridge with no groundwater recharge. The downslope end is considered to terminate in a valley with a 2 m deep water table. This conceptual model is shown schematically in figure 2.

This hillslope is represented using a 45 m deep, 1 km wide 2D domain, shown schematically in figure 4. The mesh is discretized into 100 10 m wide vertical columns and has an inclination slope of 10%. Each vertical column is discretized into 30 cells with thicknesses varying from 1 cm at the ground surface to 12 m at the bottom of the mesh, gradually increasing in size toward the lower boundary. We use finer resolution grid cells near the ground surface to capture the annual freeze/thaw dynamics in the active layer. To aid in the simulation spinup process (described below in section 2.4), a 1D column mesh was also used. The 1D mesh has the same vertical discretization as the 2D mesh.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** An exaggerated schematic representation of the hillslope domain including soil texture and boundary conditions setup.
Download figure:
Standard image High-resolution image

We differentiate between soil textures/types by applying different material properties to the cells of the mesh. The model is comprised of three layers: organic soil, mineral soil, and cobble-rich mineral soil (figure 4). Based on general observations at the Teller site, the organic soil layer is thickest at the top and bottom of the hillslope in areas dominated by tundra and grass/sedge vegetation, and thinnest in the middle of the hillslope, where shrubs predominate. In the model, the OLT is chosen to vary from 50 cm at the top of the slope to 5 cm in the middle of the slope, and 30 cm at the bottom of the slope (figure 4). The mineral soil/cobble-rich mineral soil interface is set to be constant at a 2 m depth, therefore the thickness of the mineral soil layer varies with the organic/mineral soil interface (figure 4). The organic soil thickness in the 1D mesh simulations is set to 5 cm (195 cm thick mineral soil layer), consistent with the midslope values in the 2D mesh. Thermal and hydrologic properties of the subsurface layers are summarized in table 1. As described in section 2.1, thermal conductivities in each cell are calculated at each timestep based on a formulation that considers the thermal conductivity of the phases of water and soil grains and their relative proportion present in the cell. There is no heat capacity for the bulk soil, instead we use a mixture model, which includes empirical fits for water/ice/air and a linear model for rock. The rest of the subsurface properties used in this experiment can be found on the corresponding github repository⁷ .

2.4. Experiment setup

In this section, we describe the steps required to set up the initially continuous permafrost layer. The first set of steps are performed on the 1D mesh to ease the computational burden of initialization and model spinup.

1.
Initialize hydrostatic conditions: solve Richard's equation (isothermal, above freezing) using a constant pressure of 0.54 MPa at the bottom boundary. The bottom boundary pressure was chosen to initialize the 2 m water table in the 1D column.
2.
Freeze the domain from below: initializing with pressures from step 1 and a uniform temperature of 1 °C throughout the 1D column, solve the subsurface equations with a bottom boundary temperature fixed at −3 °C and top fixed at 0 °C.
3.
Spinup the model to an annual cyclic equilibrium state: initializing with pressures and temperatures from step 2, solve the fully coupled system. The SEB is driven by the detrended meteorological data (i.e. 36 year temporal trends removed from the daily air temperature and precipitation data discussed in section 2.2). This step was repeated for 5 constant bottom boundary temperatures of 0 °C, 1 °C, 2 °C, 3 °C, and 4 °C to create different subsurface thermal regimes. These thermal regimes resulted in the corresponding permafrost thicknesses of 45, 33, 22, 20, and 17 m.

The pressure and temperatures from the 1D spinups are then assigned to each of the columns of the 2D mesh. Five 2D models are then also spun up to a cyclical thermal steady state similar to the five 1D models, i.e. using the five bottom boundary temperatures listed above. In all cases, the following boundary conditions were used: no flow at the bottom boundary, no flow along the vertical upslope boundary, hydrostatic boundary with 2 m deep water table along the vertical downslope boundary, and uniformly distributed snow precipitation (see figure 4). These simulations smooth out any temperature or pressure differences introduced by the spatially varying organic layer thickness and produce a continuous permafrost layer with spatially varying thickness in the models. For the 2D spinups, we repeat the 36 year cycle with detrended meteorological data up to three times mean annual cyclic state thermal equilibrium is achieved. Under the current setup, the sub-permafrost groundwater is still transient and slowly draining out along the vertical downslope boundary due to low permeability in the second mineral layer.

At this point, the 2D model is ready for simulations with preferentially distributed snow depth (PDSD) used for the analyses of talik formation. In these runs, we preferentially distribute snow by multiplying the snow precipitation data along the middle 240 m of the hillslope (horizontally from x = 380 m to x = 620 m, where x is the horizontal coordinate starting at the upslope end of the transect). The snow precipitation data for the rest of the hillslope remains unchanged. The multiplication factor increases linearly from 1 at x = 0 to 1.5 at x = 380 m to 3 at x = 480 m and then decreases linearly from x = 520 m back to 1.5 at x = 620 m and then to 1 at 1000 m. This formulation provides a simple representation of the effect of tall shrubs on localized increased snow accumulation.

3. Results

3.1. The effect of snow on ground temperatures

While long term air temperature trends are the dominant factor controlling ground temperature, we demonstrate the importance of snow depth on ground temperature and its ability to negate the effect of short term, seasonal air temperature in some cases. We ran the model with original (non-detrended) meteorological data to illustrate the interaction between uniform snow depth and ground temperatures. First, in figure 5, we look at the relationship between the annual average difference in temperature between the air and the ground surface (T_diff, left axis) and maximum annual snow depth (right axis). In this plot, it is apparent that although there may be some general positive correlation (Pearson correlation coefficient of 0.63) between the snow depth and temperature difference, the snow depths are much more erratic than temperature differences, and it is possible to have shallow snow depths in years with large temperature differences (e.g. 1999, 2011).

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** (A) The magenta line with dots is the maximum snow depth (SND); the blue line with triangles is the temperature difference (T_diff) calculated as a difference between mean annual ground surface (T_sur) and mean annual air (T_air) temperatures. (B) The blue line is the calculated yearly averaged subsurface temperatures (*T_g*) and the red line is the mean annual air temperature (T_air).
Download figure:
Standard image High-resolution image

Figure 5(B) shows a time series of mean annual spatially-averaged ground temperatures (T_g; left axis) and MAATs (T_air; right axis). To calculate T_g, ground temperatures are spatially averaged over the whole domain and the spatially-averaged temperatures are then averaged over each year. Figure 5(B) demonstrates that ground temperatures are not exclusively driven by the air temperature, and in some years, cold ground temperatures occur when the air temperature is relatively high (e.g. 2002 and 2013). While the first year with a low ground temperature (less than −0.2 °C in 1984) corresponds with low air temperature, the two years with the lowest ground temperatures (years 2002 and 2013) have high air temperature. Inspecting the snow depths in years directly preceding 2002 and 2013 in figure 5(A) implies that preceding years of shallow snow depth may be the cause of low ground temperatures in spite of high air temperatures. This analysis indicates that a continuous warming of the subsurface, would require continued years of sufficient snow depths to keep temperature differences high. Here, we use yearly averaged simulation results of snow depth and ground temperature. However, it is important to note that snow characteristics (Sturm et al 1997) and more importantly shoulder seasons might also affect T_diff and, therefore, ground temperatures (Gouttevin et al 2018).

3.2. Modeling of a through talik

The insulating effect of deeper snow in the middle of the hillslope drives a trend of ground warming. We investigated the formation of the talik for different permafrost thicknesses. In figure 6, we present snapshots of temperature distributions along the hillslope at the beginning (year 0), middle (year 100), and end (year 200) of the PDSD run. The first column of the plot matrix on figure 6 (year 0) shows the initial permafrost thicknesses (cool colors in the plots indicate frozen ground, while warm colors indicate thawed ground) as a result of the 2D spinup, modeled for five different temperature lower boundary conditions. Note that while the initial permafrost thickness is variable along the transects due to OLT thickness and boundary effects, we use the thickness at the center of the transects for reference.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** The snapshots of the ground temperature along the hillslope transect at different times: 0 year, 100 years, and 200 years. The 100 year and 200 year snapshots show the effect of the preferential snow depth on the ground temperatures. Each time slot has five cases representing different permafrost thicknesses simulated with detrended meteorological data.
Download figure:
Standard image High-resolution image

After 100 years, a through talik was developed only for the 17 m thick permafrost layer. The warmer temperatures beneath the talik for 20 and 22 m thick permafrost layers at 100 years indicate the warming from below that contributes to the development of a through talik at 200 years.

After 200 years, 11 m deep taliks have developed for the 45 and 33 m thick permafrost layers. The through taliks were developed for 20 and 22 m thick permafrost layers. The time required for development of a through talik varies according to the thickness of the permafrost and subsurface thermal conditions. The time lag between through talik development for 45 and 33 m thick permafrost layers was 18 years and between 20–22 m thick permafrost layers was 28 years.

The model results indicate that water movement in the sub-permafrost region increases significantly after the through taliks form. Since the upslope vertical boundary of the model is considered a hilltop without recharge, the sub-permafrost region has residual drainage at the downslope vertical boundary of around 0.3 m³ yr⁻¹, as shown at early times (80–90 years) in figure 7. After the through taliks form, the volumetric flow rates increase from 0.3–0.7 m³ yr⁻¹. This indicates that the development of a through talik connects water from the active layer with sub-permafrost water and creates new groundwater flow paths.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** The change in sub-permafrost flow after a through talik was developed. The red, blue and green curves illustrate increases in drainage for the 22, 20, and 17 m permafrost thicknesses, correspondingly. The dashed black line corresponds to the divider line in flow before and after a through talik was formed.
Download figure:
Standard image High-resolution image

To assess temporal changes in hydrothermal conditions with respect to the first year of the PDSD run, we calculate the percent difference relative to the first year according to the following formula:

$\begin{eqnarray}&&R=\displaystyle \frac{{\bar{x}}_{i}-\overline{{x}_{1}}}{\overline{{x}_{1}}}\cdot 100 \% ,\end{eqnarray} \tag{ 1 }$

where $\bar{x}$ represents the state variable of interest, and subscript 'i' and '1' represent the current year of interest and first year of the PDSD run, respectively. Here $\bar{x}$ is spatially averaged over the whole domain and then averaged over the corresponding year. Figure 8 shows spatially and annually averaged relative changes for four subsurface state variables: liquid saturation s_l, ice saturation s_i, ground temperatures T, and the total soil water mass s_m (mol) (including water in both liquid and solid states). The results of the five PDSD runs showed the highest relative increase in s_l for 45 and 33 m thick permafrost layers. The s_i shows a monotonic decrease with respect to permafrost thickness, except for 22 m permafrost layer. Relative changes in T and s_m showed overall increases.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** The relative difference between four spatially and annually averaged soil parameters at the beginning and end of the distributed snow depth run. Each bar is associated with a permafrost thickness.
Download figure:
Standard image High-resolution image

To illustrate the spatial changes within the hillslope, we plot the difference in ground temperatures (T_g) between the beginning and end of the PDSD runs, calculated for yearly averaged temperatures at each grid cell (first column in figure 9). Temperature differences indicate highest temperature increase within the areas with talik, extending in horizontal and vertical directions along the sides of the hillslope and downward beneath the talik. The cooling along the sides of 20 and 22 m thick permafrost could be due to the fact that the system for these two cases might not reach complete equilibrium during the 2D spinup.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** The first column shows the changes in ground temperatures at the beginning and end of the distributed snow run. The second column shows the relative changes in total water mass between the end and the beginning of the distributed snow run. The rows correspond to the permafrost thicknesses.
Download figure:
Standard image High-resolution image

The relative changes in s_m between the beginning and end of the PDSD runs similarly were calculated for each grid cell within the hillslope (second column in figure 9). The relative changes in s_m indicate changes in the subsurface water storage. The 45 m case shows a relative increase of s_m in the talik area and no changes in the rest of the system excluding the active layer, since the rest of the subsurface remained frozen. All other cases show increases in s_m within the talik, which correlates with the temperature changes on the first column in figure 9. Note that water continues to move slowly within the sub-permafrost layer due to residual drainage described above and shown on figure 7. The orange colored areas along the hillslope correspond to the areas that lose water, where dark violet areas correspond to the areas that gain water. Changes in s_m show groundwater movement from the top of the hillslope to the bottom. The results of the spatial changes over time shown on figure 9 were summarized on figure 8, illustrating overall changes in the 2D hillslope.

In addition to the relative changes in s_m we calculated streamlines emanating from five locations starting within the talik before and after the through talik formed to illustrate changes in the hydrologic flow regime. Figure 10(A) shows the streamlines calculated from the velocity vector field when the through talik is still forming, not allowing groundwater to penetrate to the sub-permafrost region. Figure 10(B) shows streamlines after the through talik is formed, indicating that the through talik leads to above permafrost water flowing to the sub-permafrost.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Magenta colored streamlines for five different initial starting points within the talik. The streamlines indicate the calculated flow paths for five points starting in the talik (A) before and (B) after a through talik was formed. The colorbar corresponds to the corresponding temperature distributions.
Download figure:
Standard image High-resolution image

To illustrate the evolution of the relative changes in s_m over time we show four snapshots, each 50 years apart, two before and two after a through talik formed in a 20 m thick permafrost (figure 11). The percent change in the mass of water for each grid cell with respect to year 1 was calculated according to equation (1). Here, we restrict the range in percent difference between −2% and 2% to better highlight the changes in s_m within the sub-permafrost layer. Years 50 and 100 show positive changes in s_m within the talik. The changes in the sub-permafrost layer are associated with residual water flow due to the model setup (see section 3.2). The total percent change of water in the whole hillslope R is negative, indicating that under the current setup the system loses more water overall than it gains. R changes from negative to positive for years 150 and 200 corresponding to the existence of a through talik and indicating that the water content has increased almost everywhere in the sub-permafrost layer. The negative change in s_m along the vertical upslope boundary of the hillslope shows loss of groundwater due to the no flow boundary condition.

Figure 11. Refer to the following caption and surrounding text. — **Figure 11.** Relative percent change (R) in soil moisture within each grid cell relative to the initial year. The colorbar indicates the relative difference in total water mass in each grid cell between the time indicated in the subplot and year one.
Download figure:
Standard image High-resolution image

4. Discussion

Our numerical experiments show that preferential snow distribution can lead to talik formation, which can explain the localized absence of near-surface permafrost at the Teller site. The development of a through talik depends on the permafrost thickness and localized snow accumulation to provide enough insulation (i.e. temperature difference between air and ground surface) to keep the ground warm. It is important to note that soil thermal properties, soil water content, and organic layer thickness are important parameters that can slow down or speed up the time necessary to develop a through talik. In our simulations, in addition to the boundary conditions along the bottom of the domain and mesh resolution, soil hydrothermal conditions can affect the time of the talik development. The actual timing of the talik formation in this study is conditional on our numerical setup, and future work is needed to include a sensitivity analysis to generalize these results across different ecosystems or to apply this setup to a specific geographical and landscape conditions. Development of a through talik results in a new hydrologic flow pathway that leads to interaction of surface, active layer, and sub-permafrost waters. The evolving flow paths resulted in increased sub-permafrost water storage and contributed to an increase in sub-permafrost flow. Current results improve our understanding of the influence of talik formation and evolution on subsurface hydrothermal conditions, emphasizing that the development of taliks could lead to an increase in subsurface water storage.

Snow-shrub interaction experiment near Council, Alaska suggests that areas with taller shrub vegetation collect and retain significantly more snow than the surrounding short statured tundra (Sturm et al 2005). Another snow and vegetation survey near Happy Valley site in arctic Alaska, records deepest snow under tallest and densest shrubs (Sturm et al 2001). Our research indicates that this, in turn, may lead to talik formation. O'Neill and Burn (2017) examined multi-decadal thermal effect of snow fences in Western Arctic Canada, where MAAT is −7 °C. They observed higher soil moisture and talik development near snow fences. Their numerical simulations suggest a talik developed 25 years after the fence was constructed, with a thaw depth of about 3 m, suggesting that warming will continue and might take up to 200 years to reach a new equilibrium. Our results show that even after 200 years, the simulations did not reach an equilibrium.

The present version of the model does not include vegetation; instead we control snow accumulation by preferentially distributing snow precipitation along the hillslope. While vegetation will remove some of the groundwater and alter the soil moisture distribution due to transpiration, we expect the effect of vegetation to minimally impact these results due to the abundance of water and relatively small biomass of vascular plants in these ecosystems. Moreover, near-surface thawing might trigger soil subsidence (Hinkel and Hurd 2006). The magnitude of the soil subsidence depends on the amount of near-surface ground ice, soil overburden thickness, and soil organic content. There is a wide range of literature reporting soil subsidence in the tundra (Jones et al 2013, Lui et al 2014, Schaefer et al 2015, Liljedahl et al 2016, Streletskiy et al 2017). At the Teller site, tall shrubs are located in well drained areas with thin overburden (1–2 m) underlain by competent bedrock, and may not be very susceptible to thermokarst processes. However, if subsidence occurs this would lead to constant surface water pooling, which might further accelerate permafrost degradation (i.e. talik development).

The effect of ground warming as a result of tundra shrubification has been discussed in previous studies (Lynch et al 1999, Chapin et al 2000, Cable et al 2016). However, the ramification of snow-shrub interactions is not linear and depends on multiple factors (Sturm et al 2005). In addition to the previous studies, we show that snow-shrub induced ground warming might introduce changes to the subsurface hydrology, where through taliks can provide additional flow pathways for transport of nutrients and soil carbon, increasing yearly groundwater discharge. Increases in groundwater discharge can possibly contribute to the observed increases in Arctic river discharge during the winter (Peterson et al 2002, Walvoord and Striegl 2007). However, further research is needed to better address how and to what extent the effect of through taliks can be scaled up to the watershed and larger scales.

5. Conclusions

There are two major mechanisms for natural, preferential snow distribution in the Low Arctic tundra: wind-topography interactions and interception of snow by vegetation. In this study, we omit the effect of local small-scale topography by implementing a flat, inclined surface along the hillslope within a watershed. Assuming that tall shrubs exist on the hillslope and trap more snow than surrounding tundra, deeper snow is likely to occur in that area. To better understand permafrost heterogeneity along the Arctic hillslope, we model talik development for different permafrost thermal regimes corresponding to different permafrost thicknesses, ranging from 17–45 m. Our modeling results confirm that preferentially distributed snow leads to permafrost thawing, suggesting that ground under tall and dense shrubs might either have deep talik or no permafrost. Transition of the landscape with continuous permafrost to permafrost with through taliks (i.e. discontinuous permafrost) might take from years to decades and even centuries depending on the permafrost thickness and hydrothermal conditions of the ground. Our modeling study shows that existing along the Arctic hillslope shrubs that constantly trap more snow than surrounding tundra lead to permafrost degradation even in the absence of warming. New pathways developed by a through talik allow linkages between surface water and sub-permafrost groundwater, affect changes in permafrost hydrology and lead to increasing sub-permafrost flow.

Acknowledgments

This work is part of the Next Generation Ecosystem Experiments (NGEE Arctic) project which is supported by the Office of Biological and Environmental Research in the DOE Office of Science.

Dates

Peer review information