RESEARCH/REVIEW ARTICLE

Air warming trends linked to permafrost warming in the sub-Arctic catchment of Tarfala, Sweden

Romain Pannetier^1,2 & Andrew Frampton^1,2

¹ Department of Physical Geography, Stockholm University, Svante Arrhenius väg 8, SE-106 91 Stockholm, Sweden
² Bolin Centre for Climate Research, Stockholm University, Svante Arrhenius väg 8, SE-106 91 Stockholm, Sweden

Abstract

Recent ground temperature records from the 100-m-deep borehole near the Tarfala Research Station in northern Sweden reveal that permafrost is warming at a pace consistent with the rate of measured air temperature increase at the site. Here we investigate whether air temperature increase is the main driver of the observed change in the permafrost thermal regime using a non-isothermal hydrogeological numerical model for partially frozen ground. The local site is investigated with different ground surface temperature scenarios representing different integrated effects of surficial heat attenuation processes. Results indicate that despite a short-term sensitivity to heat attenuation processes including snow conditions, the main driver of change in the permafrost thermal regime during the past decade is warming air temperatures. Additionally, the approach used here is shown to be particularly pertinent for modelling warming trends, despite limited prior knowledge of site-specific conditions and geological properties. Understanding the main driving mechanisms of changing permafrost is useful for assessing the suitability of borehole temperature records as proxies for past environmental conditions as well as for modelling possible future climatic impacts.

Keywords
Permafrost; warming ground; PACE borehole; numerical modelling; climate change.

Permafrost covers approximately 17% of the Northern Hemisphere (Zhang et al. 2003) and therefore constitutes an important terrestrial subsystem. The Global Climate Observing System considers permafrost borehole temperatures and active layer thickness as two “essential climate variables” (Bojinski et al. 2014) for monitoring climate change. Since climate is anticipated to undergo continued global warming (Solomon et al. 2007), and since the pace of this warming is faster in circumpolar regions (Hansen et al. 1999; Turner et al. 2007), permafrost is expected to undergo significant changes in the near future in terms of its distribution, thermal state and ALD.

The warming of permafrost is strongly correlated to increasing atmospheric temperatures, especially when considering longer timescales (Romanovsky et al. 2007). Yet several other factors are known to affect permafrost temperature, including incoming radiation, the aspect and slope of the ground surface, the thickness, density and distribution of snow, the vegetation cover and soil moisture (Guymon 1976; Hillel 1980; Williams & Smith 1991; Etzelmüller 2013).

In an early effort to simulate permafrost temperatures using a numerical model of heat diffusion accounting for latent heat effects of ice–liquid phase change, Kane et al. (1991) represented climatic changes by adding a linear increase to GST. Osterkamp (1984) also used increasing GST to predict a slow reduction of the permafrost thickness at Prudhoe Bay, Alaska. However, in those studies, no assumption was made regarding the driving mechanisms of the GST warming. Using a model for integrated snow cover and permafrost, that is, an extension of the SNOWPACK model, Luetschg et al. (2008) showed that the GST is more sensitive to the depth and duration of the snow cover than to the mean annual air temperature. More specifically, they claimed that the GST decreases exponentially in response to a thinner snow cover, whereas the response to a lower mean annual air temperature is merely linear.

Several studies have highlighted the complex effects of changes in winter precipitation on ground temperatures (Osterkamp & Romanovsky 1999; Stieglitz et al. 2003; Marmy et al. 2013). The accumulation of a snow cover is well known to insulate the ground during the cold season, resulting in less heat losses (Goodrich 1982). Field experiments have shown how a sudden increase in snow depth quickly results in important ground warning (Zhang et al. 1996; Zhang 2005; Hinkel & Hurd 2006). More generally, precipitation impacts ground thermal conditions. For example, Engelhardt et al. (2010) showed how increased autumn precipitation can lead to either a significant increase or decrease of the ALD depending on whether it falls in liquid or solid form. Precipitation also induces infiltration that results in heat advection (Scherler et al. 2010) and influences soil moisture content, effective thermal conductivity and storage of latent heat (Woo 2012).

Many modelling efforts have striven for an accurate representation of surface process dynamics. Hinzman et al. (1998) predicted a transient active layer thickness over a summer on a regional scale using a surface energy balance coupled with a one-dimensional subsurface model for thermal diffusion. That model was driven by regional interpolations of air temperatures, relative humidity, wind speed, incoming shortwave radiation and net radiation. Ling & Zhang (2004) presented a one-dimensional heat transfer model with a surface energy balance accounting for snow cover, which additionally required snow depth information. More recently, Endrizzi et al. (2014) developed a model of snow and soil that considers both water and heat transport together with a distributed model for snow cover, a soil freeze–thaw algorithm, as well as accounting for a subsurface water flow. That model also allows for a complex topography, which influences the energy balance by modulating the solar radiation with the incoming angle of solar radiation.

The question of surface process representation is also of interest on a broader scale. Koven et al. (2013) observed that differences between global climate model predictions of permafrost temperature and thaw were mainly due to the different representations of land-surface dynamics. Hence efforts have been made in the direction of integrating surface processes in global land-surface models. Chadburn et al. (2015) introduced a land-surface scheme that represents snow cover with multiple layers, a transient evaluation of moisture content based on the Richards equation and Brooks and Corey retention curves, as well as a representation of the organic soil properties. These physically based models allow comprehensive representation of the surface processes and their interactions, so that the system’s sensitivity to different changes can be explored and accurate predictions made about the upcoming changes in permafrost, assuming the availability of reliable inputs. These are also very data-intensive schemes; Atchley et al. (2015) give an account of the complex procedure involving recursive calibration and model structural adjustment that may be necessary to parametrize a thermal hydrology model with representation of surface processes (Coon et al. 2013). Application of physically based models can therefore be limited by lack of data.

A common empirical method uses so-called “n-factors” to determine GST from air temperature in a considerably more data-conservative manner. It substitutes energy balance calculations with empirical approximations based on air and ground temperature observations. As such, the factors account for the near-surface heat attenuation processes of climate dynamics in a lumped way (Lunardini 1978). Typically, two n-factors are used to represent the air–ground relationship at a given site, one for freezing conditions and one for thawing conditions (Riseborough et al. 2008). The freezing n-factor n_f is defined as the ratio between the freezing degree-day sum of the soil surface and the freezing degree-day sum of the air, where the freezing degree-day sum is defined as the integrated departure from 0°C over days with negative temperatures. Thawing n-factors n_t are defined similarly for positive temperatures (Lunardini 1978; Smith & Riseborough 1996; Riseborough et al. 2008; Farbrot et al. 2013).

Juliussen & Humlum (2007) showed that the effect of solar radiations and snow depth could be successfully characterized by n_t and n_f, respectively. Yet n-factors do not discriminate between the different physical processes that control the relationship between air temperatures and GST (Karunaratne & Burn 2003). Furthermore, since snow depth exerts a strong control on n_f, and since the spatial and temporal variation of snow conditions is high in sub-Arctic mountainous environments, n-factors can have a strong spatial and temporal variability (Karunaratne & Burn 2003; Etzelmüller et al. 2007; Juliussen & Humlum 2007). An inherent limitation related to the empirical nature of n-factors is therefore that they are not necessarily suitable for extrapolation to other sites, nor to past or future climates.

Following the recent emergence of models for coupled groundwater and heat flow for permafrost environments (McKenzie et al. 2007; Painter 2011; Hammond et al. 2012), several studies have focused on climate warming effects on groundwater flow dynamics in permafrost systems at decade to century timescale. Ge et al. (2011) used air temperature to drive their model and to simulate anticipated climate changes, while Bense et al. (2009), Frampton et al. (2011), Bense et al. (2012), Frampton et al. (2013) and Frampton & Destouni (2015) performed simulations driven with GST whose increases are consistent with predictions of air temperature warming.

Jonsell et al. (2013) documented substantial changes in the thermal regime of permafrost at the Tarfala site, northern Sweden, and calculated linear warming trends at different depths at and below the depth of zero annual amplitude, based on borehole temperature measurements between 2001 and 2011. They noted a consistency between increases in air and ground temperatures measured in the Tarfala PACE borehole over this time period. Considering the conditions and surroundings of this borehole, specifically its location on a ridge, groundwater flow is expected to be minimal and snow cover is assumed to be thin (Isaksen et al. 2001; Isaksen et al. 2007; Christiansen et al. 2010). In this context, Jonsell et al. (2013) highlighted the question of whether the changes in ground thermal regime could solely be explained by changes in air temperature. Regression analysis suggests that thermal conditions in the ground are correlated with air temperatures (Christiansen et al. 2010). If this indeed is the case, it could provide further support for the suitability of the Tarfala site to be used in the context of retrieving palaeoclimatic temperature records, which is one of its stated objectives as part of the PACE programme.

The aim of this study is to investigate the hypothesis that air temperature increase is the main cause of recent permafrost warming at the location of the PACE borehole near Tarfala, northern Sweden. An additional aim is to identify the potential influence surficial heat attenuation processes may have on the development of subsurface temperatures, with a particular focus on trends and the changes in the permafrost thermal regime. A numerical model which solves for the system of coupled differential equations governing heat and moisture balances in partially frozen ground is applied, and then the resulting modelled ground temperatures are compared against measured borehole temperatures. Three scenarios are considered, in which surficial heat attenuation processes are either neglected, accounted for implicitly, or approximated and maintained constant over the study period, which serves to enable distinction between the respective influences of air temperature and surficial heat attenuation processes on ground temperature changes.

Method and site data

Site description

The Tarfala Research Station (67°54.7 N, 18°36.7 E) is located at 1135 m a.s.l. in a sub-Arctic, mountainous environment (Fig. 1). A long-term monitoring programme has been conducted at this site since 1965, primarily with focus on the mass balance of the glaciers, concurrently with meteorological measurements to allow studies of glacier–climate interactions. A description of the site and overview of available data are provided by Jonsell et al. (2013); in the following section, a summary of the data relevant for this study is provided.

Fig. 1 Overview map of the Tarfala valley, indicating the location of the Tarfala Research Station (red star), the position of gauges along the Tarfala stream (blue squares) and the three PACE boreholes (pentagons). The topographical catchment boundary is marked by the red line. Elevation difference between isolines is 50 m.

The site hosts a total of three boreholes drilled as part of the PACE project. The main borehole, denoted as PACE-1 (Fig. 1), is located on the saddle (67°55′ N, 18°38′ E), called Tarfalaryggen, close to the topographic water divide of the Tarfala catchment at elevation 1551 m a.s.l. It was drilled in 1999 down to a depth of about 100 m, has a diameter of 16 cm and is equipped with casing and a chain of 44 006 model thermistors (Yellow Springs Instruments; 0.02°C relative accuracy) whose depths are reported in Table 1 (Harris et al. 2001). At this location, the seasonal weather patterns exert an important influence on ground temperatures down to the depth of zero annual amplitude at ca. 20 m. Here, zero amplitude is defined as an amplitude lower than 0.1, consistently with the definition of Isaksen et al. (2001). The zone beyond this threshold is henceforth referred to as the “geothermal zone” (Kurylyk et al. 2014).

A 15-m-deep control borehole, denoted as PACE-2 (Fig. 1) was drilled near PACE-1. A record of the temperature time series along the profile is available since March 2000 (Isaksen et al. 2001). This data set, for the time period 2001–05, has later been described and analysed by Isaksen et al. (2007). A third borehole (PACE-3) was drilled in 2008 in the valley close to the Tarfala Research Station (Fig. 1). It was drilled down to a depth of 15 m and installed with 18 thermistors, but because of measurement complications contains only a partial record of ground temperatures.

The PACE-1 drill-site is barren, with almost no vegetation except for scarce lichen, mosses and few occurrences of the vascular plant Ranunculus glacialis. It is covered with an unconsolidated regolith that shows on the surface and that is affected by sorted polygons (Isaksen et al. 2001). The regolith has been reported to be 4 m deep and the underlying bedrock, which consists of amphibolite (Andréasson & Gee 1989) and is fractured down to a depth of about 10–15 m (Isaksen et al. 2001).

Daily air temperature averages measured in the valley, that is, near the Research Station, are available since 1965 and manual measurements of summer precipitation have been conducted since 1980. However, the precipitation record is limited to a varying time period for which the station is manned, and winter precipitation is only recorded as snow depth over the monitored glaciers. This introduces difficulties in assessing both seasonal and total annual precipitation rates.

Data processing

Air temperature at the PACE-1 borehole location is measured at a height of 2 m above ground and with a six-hour resolution. Measurements have been conducted since the year 2000; however, because of technical problems this record is subject to several periods without measurements. Therefore, in this study, data gaps are filled by deriving temperatures from measurements obtained from the air temperature measurements at the Tarfala station. This is done by using the lapse rate calculated by Jonsell et al. (2013), which gives air temperatures at PACE-1 as a function of air temperatures at the Tarfala Research Station. The same lapse rate is also used to extrapolate air temperatures at the PACE-1 location for time periods prior to installation of the air temperature sensor. In order to differentiate between the impact of the warming air temperatures and the influence of heat attenuation processes acting on the ground surface, the numerical model is run with three different temperature series as surface boundary condition. One model scenario is based on directly applying the air temperature record, and is denoted as T-Air. Another model scenario is based the GST record, denoted as T-Gnd, which corresponds to the most shallow thermistor (0.2 m) in the borehole. T-Gnd therefore reflects the combined influences of air temperature and transient heat attenuation processes occurring on and in the ground surface down to 0.2 m depth. The third model scenario and boundary condition, denoted as T-Est, consists of estimated ground temperatures which are obtained by combining air temperature series with two linear transfer functions. By keeping these functions constant over time, the T-Est model case has a boundary condition whose changes only follow the evolution of the surface air temperature. Hence, for this purpose, an empirical function is more appropriate than a process-based energy-balance model. It is also convenient since such a scheme is obtainable even in the absence of snow depth data. The estimated ground temperature T-Est is obtained from the air temperature time series T-Air(t) by:

for summer t_s and winter times t_w, respectively, and where the coefficients of the linear fits are obtained from field data (Fig. 2). These linear equations correspond to a fully empirical method for relating measured air temperatures to an estimate of corresponding ground temperatures applicable to the specific site investigated. It is adopted because it maintains this relationship constant over time, corresponding the time frame of available data (2001–11). Similar linear fits between air and GST measurements have previously been used by Zhang et al. (1996) and Zhang et al. (1997) for snow free conditions. There is also a strong similarity with methods based on n-factors (Lunardini 1978), which however typically integrate temperatures over time, distinguishing between negative (freezing) versus positive (thawing) temperatures, thereby resulting in varying durations of freeze–thaw periods as changes in the temperature record occur.

Fig. 2 Daily observations of ground temperatures (at 0.2 m depth) plotted against corresponding air temperatures at PACE-1 over the period 2000–2011. The data points are sorted in two groups, summer (green) and winter (blue), from which two linear best-fits (red lines) are obtained, and used as estimates of surface ground temperatures.

Furthermore, to smooth variability in temperature time series data the following 11-day central moving average algorithm is used:

where 〈T〉_i is the resulting averaged temperature for day i, T i is the daily temperatures and N is the total number of days in the temperature record. This smoothing is applied to both the air and ground temperature time series, as well as to the estimated ground temperature representations T-Est.

Simulation configurations

The numerical model represents a one-dimensional column of the subsurface ground in the vicinity of the PACE-1 borehole. The top of the model, z=0, coincides with the ground surface at 1551 m a.s.l. and the bottom extends 100 m below the bottom of the borehole, to a depth of −200 m.

The top of the domain is assigned a constant liquid pressure corresponding to atmospheric pressure, implying zero capillary pressure at the ground surface. This maintains a steady water table at the top of the domain, and fully saturated conditions within the domain, such that only liquid water and/or ice can occupy the pore space. The bottom of the domain is assigned a no-flow boundary for water and set to a constant temperature T=–2.1°C corresponding to a linear extrapolation of temperatures measured at the bottom of the PACE-1 borehole to 200 m depth. Lateral sides are assigned no-flow boundaries for both water and heat.

Surface temperature is represented by a transient temperature time series obtained from data, based on an 11-day moving average, Eqn. 2 and imposed with a daily temporal resolution. The following three cases for surface temperature representation are considered and form the basis of the comparisons conducted in this study.

In the first case, denoted as T-Air, surface temperatures are directly assigned the 11-day moving average of air temperature time series. This scenario implies an assumption of air temperature coinciding with GST and hence does not account for any surficial processes attenuating heat transfer to the ground.

In the second case, T-Gnd, surface temperatures are assigned the 11-day moving average of ground temperatures obtained from the PACE-1 borehole. Here the top-most sensor is used, which is located at a depth of 0.2 m below the ground surface. Hence, this scenario serves to implicitly account for the net effects surficial processes can have on heat transfer, for example, snow insulation, radiation and evapotranspiration, as well as heat transfer effects occurring in the upper 0.2 m ground layer.

In the third case, T-Est, surface temperatures are obtained by converting air temperature to an estimated ground temperature using linear correlations obtained by Eqn. 1. This scenario assumes that the net effect of surficial heat attenuation processes will remain unchanged over the studied period. The T-Est boundary condition is derived from air temperatures, using two empirical linear functions, which represent the average thermal dampening effect of the surface. The assumption made here is that if this scenario can successfully reproduce the observed temperature trends in the permafrost, then it indicates that this evolution is primarily induced by air temperature.

In order to initialize the computational domain, suitable initial conditions which are physically consistent with the imposed boundary conditions need to be assigned. This involves obtaining liquid pressure and temperature fields and associated ice–liquid phase states in the subsurface. The necessary procedure to achieve this is divided into a sequence of four steps consistent with the recommendations provided by Karra et al. (2014). First, an unfrozen hydrostatic system with water table at z=0m is defined. Then a cryostatic system is obtained by cooling from below with the assigned constant bottom temperature T=–2.1°C until the entire system becomes fully frozen and stabilizes at this temperature. During this stage, small quantities of liquid water may flow out of the domain through the top surface as the system undergoes the liquid-to-ice phase change and associated volume expansion and liquid pressure changes.

The third step is to impose an annually periodic temperature signal representative of seasonal variability, as well as to establish a climatic history of the temperature variation with depth. To achieve the former, a sinusoidal representation of the earliest available annually-complete air temperature time series, namely from year 1965, is run iteratively. A periodic steady-state is assumed to be reached when the temperature difference in the domain between two consecutive iterations reaches below a threshold value of 0.01°C.

The fourth and last part of the initialization consist of running a simulation for 35 years using a temperature time series representative of the period 1966–2001. The boundary conditions used for this are based on air temperatures for the T-Air case, and estimated GST, using the relationships in Eqn. 1, for the cases T-Gnd and T-Est. This initial run allows integrating the influence of the recent climate history on the initial temperature profile.

The four steps described above constitute the model spin-up. The simulation itself is run over the period 2001–11, for which reliable PACE-1 borehole measurements are available. Over this period, comparison between the simulated and measured thermal state of the entire borehole.

Material properties

The model consists of three layers with material properties representing different geological strata consistent with observations at the Tarfala site. The uppermost layer consists of regolith (0–4 m), followed by fractured (4–14 m) and solid bedrock (14 m and below). Porosity, density and permeability values for each layer are obtained from the literature (Bear 1979; Dingman 2002), where regolith and amphibolite are assigned values corresponding to coarse sand and crystalline rocks, respectively (Table 2).

Effective thermal diffusivity D of the sound bedrock is estimated from site data using the method of amplitude attenuation (Vonder Mühll & Haeberli 1990; Isaksen et al. 2000), described by:

where p is the period of the signal (one year, in the present case) and m is the slope of a linear fit obtained by plotting the natural logarithm of the annual amplitude ln(A(z)) against depth z (Fig. 3). The calculation has been carried out at depths of 10–20 m for year 2009. Within this depth range, the bedrock can be regarded as relatively homogeneous and a pattern of seasonal temperature variation can be observed. Also, at those depths, the temperature time series are relatively smooth such that the differences between annual maximum and minimum temperatures are used as a sufficiently accurate estimate of amplitude. The evaluation of diffusivity has been confirmed with an accuracy of three significant digits when repeated for 2006.

Fig. 3 (a) The attenuation of the annual amplitude of ground temperatures. The amplitude A is defined as the difference between the maximum and minimum temperatures observed in one period (one year) at a given depth z. Annual minima and maxima are emphasized by the horizontal dashed and dashed-dotted lines, respectively. (b) The slope m obtained by fitting a line to the natural logarithm of the annual amplitude ln(A(z)) against depth z. The diffusivity D(m²s⁻¹) is then calculated as D=π·m²/p where p is the period.

A generic value for thermal conductivity has been assigned to the amphibolite bedrock: κ_bedrock=2.897 W·K⁻¹·m⁻¹ (Robertson 1988). The heat capacity of the bedrock is then derived by c_bedrock=κ_bedrock/D_bedrock, where the value obtained is C_bedrock=528 J·kg⁻¹·K⁻¹.

To specify the thermal properties of a soil in the numerical model PFLOTRAN-ICE (for model description, see subsequent section below), four parameters need to be specified: the heat capacity of the soil constituent C_grain and the bulk conductivity of the thermal conductivity of the soil under dry (κ_dry), fully water-saturated (κ_wet,u) and fully ice-saturated conditions (κ_wet,f). The latter three quantities are used to determine the effective thermal conductivity using Eqn. 5. Since the soil is a regolith, its grain components are similar in nature to the underlying bedrock, hence the value of C_grain is assimilated to C_bedrock and identically applied to all three layers.

Observe that whereas C_bedrock is a property of the grain material, κ_dry, κ_wet,u and κ_wet,f are properties of the bulk soil. For those properties, porosity and the nature of the pore-filling constituent come into consideration. For each layer, k_dry, k_wet,u and k_wet,f are calculated based on the assigned porosity and standard thermal conductivity values for air, liquid water and ice (Table 2) in combination with the following equation that relates the bulk thermal conductivity of soil κ_bulk to the conductivities κ_i and volumetric fractions f(i) of its components i (Woo 2012):

This equation accounts for porosity, since porosity affects the respective volumetric fractions f(i) of the bulk soil occupied by the rock, air, water and ice components.

Model overview

The numerical simulation experiments are conducted using the well-established code PFLOTRAN (Hammond et al. 2012), which solves multiphase flow, heat and transport equations and has had a wide range of applications within hydrogeology (Gardner et al. 2015; Tutolo et al. 2015). In this study, the recently developed extension PFLOTRAN-ICE for problems involving heat and water flow in partially frozen ground is used. Details of the PFLOTRAN-ICE model, including governing equations, constitutive relations and the numerical solution approach used, are described by Karra et al. (2014); in the following only the aspects of the model which are relevant to this study are presented.

While PFLOTRAN-ICE can simulate heat diffusion, groundwater flow and the ensuing advection of heat, considering three phases of water, the potential complexity of the model is limited by design in the present study. The surface boundary condition is set to a constant value of one atmosphere to ensure that the system remains saturated with liquid water and/or ice. Further, since a one-dimensional vertical domain is considered, no lateral water flow and, hence, no heat advection takes place. The dominant physical process considered in our study is therefore heat conduction with latent heat transfer in a porous media partially saturated with ice and/or liquid water, and here the representations for effective thermal conductivity and phase partitioning play a significant role.

Phase partitioning in the porous media and latent heat transfer effects are directly accounted for by thermodynamic constraints derived from the Clausius-Clapeyron relation (Painter & Karra 2014). Ice is considered immobile and reduces the available pore space for liquid water and vapour. This is combined with retention curves, based on the van Genuchten model (van Genuchten 1980) and generalized to give the respective saturation of ice, liquid water and vapour (if present) in the pore space.

Effective thermal conductivity κ has a functional dependence on the phase partitioning in the pore space, expressed as:

where κ_wet,f, κ_wet,u and κ_dry are model inputs that are calculated during preprocessing using Eqn. 4. κ_wet,f, κ_wet,u are ice-saturated (frozen) and liquid-saturated (unfrozen) soil thermal conductivities, respectively (Painter 2011). κ_dry is the thermal conductivity of the soil under dry conditions. Ke_u and Ke_f are so-called Kersten numbers, that is, the ratios of partially to fully saturated thermal conductivity for unfrozen and frozen conditions (Andersland & Ladanyi 1994). In PFLOTRAN-ICE, these are implemented as power functions of the phase saturations Ke_u=(S_l)^αu and Ke_f=(S_i)^αf, where S_l and S_i correspond to the liquid and ice saturation, respectively. The α_u and α_f are unfrozen and frozen exponent-fractions, respectively; these are set to α_u=0.45 and α_f=0.95, consistently with the values used by Karra et al. (2014).

Note that in general no model calibration is performed. Most values assigned for the various model parameters (see Table 2) are generic and/or literature values or, when possible, inferred from available field data (such as thermal diffusivity D). A sensitivity study investigating effects of realistic ranges of heat capacity and porosity values is shown in the supplementary material.

Results and discussion

The thermal state of the subsurface

The simulated ground temperatures and measured borehole temperatures are presented in Figs. 4 and 5. Simulated ground temperatures generally follow the measured temperature fluctuations, with the limitation that all simulation cases tend to underestimate measured temperatures.

Fig. 4 Temperature time series at 2, 5, 10 and 20 m depth for the time period 2000 to 2011, as measured by PACE-1 borehole thermistors (solid black), and compared with corresponding temperatures obtained from simulation cases T-Air (dashed-dot blue), T-Est (dotted green) and T-Gnd (dashed red).

Fig. 5 Measured and simulated temperatures in the shallow ground. Each column of pixels in the colour-charts represents the temperature profile for one day. (a) Observed ground temperatures T_obs(t,z). Three consecutive panes are used to present each scenario: (b), (e) and (h) show the temperature boundary conditions, for cases T-Air (blue), T-Est (green) and T-Gnd (red), with the black thin line representing the daily ground temperature from the 0.2 m thermistor. Panes (c), (f) and (i) are the simulated temperatures T_sim(t,z) and panes (d), (g) and (j) are differences between observations and simulations ΔT(t,z)=T_sim(t,z)−T_obs(t,z). The bottom pane (k) shows the depth integrated absolute as a function of time.

From the time series (Fig. 4), it is apparent that T-Air yields the largest underestimation, with temperatures being too cold at all depths. Comparative examination of the series at different depths reveals that cold winter temperatures are responsible for this discrepancy. Indeed, at 0.2 m depth (Fig. 4a) all three simulations cases show similar temperatures during summer, while the T-Air case departs notably from the others during winters. Those cold temperatures eventually influence the entire year, as seasonal fluctuations even out with depth. This indicates that the snow cover at PACE-1 has a non-negligible warming effect on the ground. Considering the site’s thin snow cover, an increase in albedo and emissivity could have counteracted the insulation (Zhang 2005). Indeed, those enhanced properties result in a cooling at the surface of snow, which can cool the ground if it is not compensated by the insulating effect of the snow cover (Zhang 2005).

Generally, the cases T-Gnd and T-Est perform better. At 5 m depth (Fig. 4), the major differences with observations occur during winter. The T-Gnd case shows larger interannual fluctuations than the other cases, namely the T-Gnd time series can sometimes be warmer and sometimes colder than the T-Est case. During the cold season, T-Gnd yields colder temperatures than T-Est at 5 and 10 m for the years 2002, 2007 and 2008, and warmer for the years 2004 and 2006. However, as depth increases, both simulations behave increasingly similarly. At 20 m depth, for instance, they appear as essentially the same, consistently about 0.6°C colder than observations. This difference further decreases to about 0.3°C at 100 m depth (not shown). The annual amplitudes of the simulated temperature signals are slightly larger than those for observations, and the simulated temperatures are almost in-phase with observations, only preceding the observations by a couple of weeks.

A detailed view of the thermal dynamics occurring in the shallow ground, with focus on depths ranging between 0 and 2 m, is shown for the time period 2005–10 (Fig. 5). The first pane of each series of three shows the imposed surface temperature boundary condition for that simulation case. The second panes show simulated temperatures T_sim(t,z), and the third panes show the differences between observations and simulations defined as ΔT(t,z)=T_sim(t,z)−T_obs(t,z). The colours are to be interpreted as ranging from blue, for model underestimation of temperatures, to red, where temperatures are overestimated. It appears that underestimation is already the general tendency in the upper layers of the ground. For cases T-Air and T-Gnd, with few exceptions, ground temperatures are systematically underestimated below 1.0 m depth. This suggests that the main limitation of the present modelling effort resides in the parametrization and/or the process representation of the active layer. For case T-Gnd, the underestimation is greatest when thaw conditions prevail. Underestimations can also be observed over the upper 1.5 m of soil at the onset of the thaw season, when the surface boundary temperature passes 0°C for the first time. Although they are intense, those discrepancies are short lived. Marked overestimations occur every year in the consecutive period, which is characterized by refreezing. Higher than observed temperatures during autumn refreezing are pervasive in all simulation cases. This overestimation is succeeded by another period of underestimation once all soil water has refrozen in the simulation.

The seasonal variability of the temperature signal reduces with depth, so that for both field measurements and simulations it has reached the depth of zero seasonal amplitude at 20 m. At those depths, where the signal is almost free from seasonal fluctuations, the warming trends become apparent.

Inferring a snow-induced warming trend

Figure 6 reveals sub-decadal patterns in the temperature evolution, such as a near stagnation in measured temperatures between 2005 and 2007 and a subsequent period of fast warming from 2007 to 2009. During those years of fast warming, the temperature time series shows a rupture from the sinusoidal pattern of the seasonal cycle to adopt a more step-like appearance; that is, seasonal increases in temperatures are no longer followed by a subsequent cooling but rather by a near stagnation. This suggests that mild winter ground temperatures are responsible for this fast temperature increase.

Fig. 6 Ground temperature evolution at 20 m depth. The oscillating curves show the daily ground temperatures. Diamonds show annual averages and the linear curves are trend lines, obtained by linear fit to those averages. The shaded area represents the 95% confidence interval of the estimated trend slope.

While the 2007–09 event of fast permafrost warming is captured well by the T-Gnd simulation (Fig. 6) it is not discernible for the T-Est case. Further, the different behaviours of the simulations T-Est and T-Gnd with regard to this event is also visible in the active layer. The T-Est case yields underestimations of temperatures for the winters of 2006/07 and 2007/08 (Fig. 5g); taking into account the approximate 11-month lag required for a temperature signal at the surface to reach 20 m depth, those are the winters that coincide with the fast temperature increase occurring in the geothermal zone between 2007 and 2009. However, for the years immediately before and after this event (i.e., winters 2005/06 and 2008/09), winter temperatures are largely overestimated. This compares to the T-Gnd case (Fig. 5j), where winter temperatures are consistently underestimated during the entire simulation period.

The fact that T-Est is unable to render the 2007–09 warming event highlights the fact that no fast air temperature warming is occurring for those years, which otherwise would be reflected in the T-Est boundary conditions. Rather it must come from a change in the surface processes, whose influences only affect the T-Gnd boundary conditions. Since the temperature increase appears to originate from mild winter GST, and considering the magnitude of those changes, it is reasonable to assume that this short event of quick ground temperature increase is the result of two consecutive years with a larger than usual snow cover. This interpretation is consistent with the fact that the North Atlantic Oscillation index was high in those years, considering that a high North Atlantic Oscillation is usually associated with mild and wet winters in Scandinavia.

The presence of a period of fast temperature increase characterized by mild winters (2007–09) highlights the important control winter precipitations exerts on the thermal regime, and how an increase in snow cover concurrent to a warming of the air could contribute to an enhanced permafrost warming. Therefore, the question of whether the ground temperature warming is exclusively an effect of the air temperature warming appears highly relevant.

Permafrost warming trends

To distinguish between direct effects of air temperature warming on ground temperatures and warming due to surface heat attenuation processes, most importantly snow cover, the case T-Gnd is compared against the case T-Est. Recall that the surface temperature boundary condition in T-Est is based on linear correlations between measured air and uppermost borehole temperatures; hence, a potential increase in snow insulation at the site would not be accounted for in the T-Est case. Had such an increase occurred, linear warming trends in the T-Est case would be flatter than trends in the surface temperatures record, that is, in case T-Gnd.

The warming experienced by the permafrost at 20 m depth is quantified by fitting a linear trend curve to the mean annual temperatures. This warming trend line is shown in Fig. 6, together with the 95% confidence interval of the estimated slope using the t-test. All resulting trend slope values and their confidence intervals are reported in Table 3. It appears that the simulations driven by surface ground temperatures, both estimated and observed, are able to reproduce the observe temperature trends; that is, their trends are not significantly different from the trends in the temperature record (Table 3). The T-Air case, in contrast, exhibits a significantly faster warming pace. Similar linear trends are calculated for all thermistor depths in the PACE-1 borehole (Fig. 7). The simulation cases T-Est and T-Gnd are able to estimate ground temperature trends reasonably well, with T-Est giving the best trends over most of the profile. The simulated trends for the T-Gnd and T-Est cases behave similarly in the geothermal zone and, at most depths, fall within the confidence interval of the observed trend.

Fig. 7 Slope of linear warming trends at different depth. The shaded area represents the 95% confidence intervals of the slopes obtained by linear best fit to measured data. The case T-Air (blue), T-Est (green) and T-Gnd (red) are represented, as well as a modification of case T-Gnd (cyan), which only differ in that it has a higher heat capacity C_r=600 J·kg⁻¹·K⁻¹.

With a slight increase in heat capacity and thermal conductivity, the simulation case T-Gnd reproduces observed trends within the confidence interval at all depths (Fig. 7), with the exception of the zone around the 95 m sensor, which displays somewhat odd behaviour. Here, heat capacity has been increased from 528 to 600J·kg⁻¹·K⁻¹, of the same proportion, so as to keep thermal diffusivity constant. The pertinence of such a parametrization is supported by the fact that an increase in heat capacity and thermal conductivity results in a reduced seasonal amplitude and a smaller lag shift. Both of these aspects of the ground temperature signal are in better agreement with measurements after this adjustment.

The performance of case T-Est shows that the observed subsurface warming trends can be reproduced reasonably well without directly accounting for changes in surficial heat attenuation processes. It appears therefore that the increase in air temperature is responsible for the changes in the permafrost thermal regime, as hypothesized by Jonsell et al. (2013). However, it was inferred in the preceding section that the subsurface temperatures are sensitive to interannual variations, together with an increase in thermal conductivity snow variations. The fact that warming trends could have been simulated without representing the snow cover dynamics leads us to conclude that no consistent trend in snow depth or, more accurately, in the insulating effects of the snow cover, has taken place at the site. This is consistent with observations by Farbrot et al. (2013), who found low interannual variability of n-factors at sites with scarce vegetation. This, however, does not necessarily exclude the possibility that changes in regional or even local precipitation has occurred. Since the borehole is located on a ridge, wind-driven snow movement is likely to play a major role for snow distribution. It is therefore likely that wind would maintain a thin snow cover, even in a context of increasing winter precipitation.

This observation that air temperatures control the permafrost thermal regime is encouraging for further modelling efforts. Indeed, if the absence of a trend in snow effects persists, the linear models used to obtain T-Est could be adopted in modelling future scenarios, subject to the assumption of stable surficial heat attenuation processes. This also supports the conclusions by Isaksen et al. (2001) that the data from the PACE-1 borehole are appropriate for palaeoclimatic analysis based on inversion modelling, which was one of the initial intentions for the use of this data (Sollid et al. 2000). This result, namely that the evolution of permafrost temperatures is mainly controlled by air temperature, can strictly speaking only be asserted for this site studied and for the period of currently available data (2000–2011).

The thermal gradient and its inflection point

Another measure of the rate of change in permafrost temperature is the evolution of the depth to the inflection point between negative and positive temperature gradients. Harris et al. (2009) have previously reported the depth of this inflection point and Jonsell et al. (2013) documented its steady lowering by 20 m between 2001 and 2011.

At the beginning of the monitoring period, in 2000, temperatures in the geothermal zone decreased with depth (negative gradient) down to 25 m below the ground surface, and increased (positive gradient) at greater depths. This point of inflection between gradients of different signs is highlighted for each year in Fig. 8. By taking year 2000 into account, this inflection point was even found to have migrated from ca. 25 m to ca. 50 m depth, that is a lowering of ca. 25 m in total.

Fig. 8 Temperature profiles on 7 June each year, temperatures are from (a) the borehole and (b) the simulation case T-Gnd.

Simulation T-Gnd presented in the lower two panes of Fig. 8 shows a similar lowering by 25 m, that is, from 27.5 m to 52.5 m. This decrease has a similar inverse exponential aspect as the observed inflection point, hence capturing a change consistent with observations and in agreement with the observed and simulated warming trends.

Active layer dynamics

Notably two aspects of the dynamics in the active layer can be discussed through comparisons between observed and simulated ground temperatures. One aspect is that measured temperatures between 0 and 2.5 m depth remain near the freezing point during the entire autumn and part of the winter, sometimes persisting as late as February (Fig. 5a). This phenomenon corresponds to latent heat transfer during phase change, resulting in temperatures remaining close to 0°C, and is henceforth referred to as “zero-curtain effect.” Here, the simulations typically underestimate the duration of the zero-curtain, leading to earlier drops of temperatures after the complete refreezing of the active layer. The importance and persistence of this effect allows us to infer a high soil water content at the site. The other aspect is that although the interannual changes in ALD are generally consistent with observations, simulated ALD are generally too shallow (Table 4).

It seems reasonable to assume both of these discrepancies occurring in the different simulations are related. On the one hand, a too shallow active layer should arguably lead to a shortened duration of the zero-curtain effect, since a reduced column of soil containing liquid water refreezes faster. On the other hand, an early onset of refreezing and subsequent drop in ground temperatures enhances cooling of the ground, which should inhibit the thawing of the active layer during the subsequent summer. Therefore, simulating a deeper active layer and hence prolonging the zero-curtain effect should bring temperature outputs closer to measured temperatures. Conversely, a general increase in temperature would logically result in additional thawing, and therefore in a deeper active layer.

The inaccuracies in ALD and the related zero-curtain effect, are therefore connected to the more general issue of temperature underestimation. Additional simulation cases, not presented here, have allowed ruling out several possible hypothesis for the underestimation of ground temperatures. These are discussed in the following paragraphs.

Uncertainties regarding the initial conditions have been investigated by increasing the number of annual iterations used for the spin-up, and by increasing the initial temperature of the system by 0.5°C. The former attempt yielded negligible changes, while in the latter case the additional heat initially presenting the system was rapidly dissipated in the upper 35 m, eventually reproducing the same underestimations as in the T-Gnd case.

Increasing the heat capacity together with thermal conductivity (thus maintaining thermal diffusivity constant) resulted in a slight reduction in cumulative differences, but at the cost of a too small annual amplitude in ground temperature, as can be seen from temperature time series at different depths (Fig. 4). Decreasing porosity also led to a general improvement of temperature estimations, but further reduced the duration and intensity of the zero-curtain effect. Further description and analysis of these two parameters are available in the supplementary material.

The temporal resolution used for the surface temperature boundary condition has also been studied as a potential source of bias. In particular, it has been conjectured that by taking diurnal temperature cycles into account, the number of freeze–thaw cycles would increase and transform the dynamics of latent heat fluxes in the active layer. However, a modified version of simulation T-Gnd, using surface ground temperature data with a six-hour resolution, failed in prolonging the zero-curtain effect and did not improve the general temperature underestimation. Ground temperature trends also remained virtually unchanged. The refined temporal resolution merely led to quicker temperature fluctuations in the near surface, down to about 1.0 m.

Since neither the different model variants nor the abovementioned parametrization attempts have notably improved the model performance, it seems appropriate to question the suitability of the current model structure for modelling the PACE-1 borehole. Perhaps the main model simplification made in the current study is the use of fully saturated conditions. Such conditions are not expected for the PACE-1 site, since it is located on a ridge, and where the upper parts of the ground consists of relatively permeable and coarse material. The presence of air in an unsaturated pore space should contribute to a reduced heat diffusion, since the thermal conductivity of air is notably lower than that of liquid water and ice (O’Donnell et al. 2009). However, the unsaturated soil should then also have a reduced effective heat capacity, that is a decreased ability to store and release heat, and the reduced total moisture content should cause a reduction in latent heat transfer. Hence the combined effect should result in a higher effective thermal diffusivity when temperature fluctuates around the freezing point.

Introducing infiltration in the model structure, may also improve the temperature underestimation observed at the onset of the thaw season. The heat advection resulting from the infiltration offspring meltwater could be at the origin of the deep-going temperature jump observed in early summer. This could in turn be an important control in the determination of the ALD. In such a model, a deeper ALD would have repercussions on the soil moisture, could influence the zero-curtain effect as discussed above, yet exerting interesting feedback on temperatures. Hence increasing the model complexity by allowing for unsaturated conditions and infiltration would allow to capture second-order processes, and the complex interplay between water and heat could explain some of the finer details of the observed ground temperature record.

However, to address this aspect additional site data are needed; in particular measurements of soil moisture and water table height in the active layer and how they changes over time, combined with suitable parametrization of retention curves for the soil layer representation of the upper region of the subsurface. Also, additional hydro-meteorological information—missing at the PACE-1 site—would be extremely useful, especially in regards to infiltration rates and variability over the year, accounting for snowmelt infiltration versus rainfall.

Conclusion

This study shows that air temperature change is the main driver for the changes in the permafrost thermal regime observed at Tarfalaryggen, Sweden. Surface temperature attenuation has remained essentially constant for the duration of the available time series, which reveals the absence of consistent changes in snow conditions at the site over the investigated time period (2001–2011). This result provides support for the assumption that the borehole’s temperatures are useful in the context of retrieving information about recent past climate conditions. It also gives confidence in the applicability of using near-future projected air temperature changes for modelling subsequent development of the subsurface permafrost regime. Finally, considering surficial heat attenuation processes has been further confirmed to be important when modelling permafrost temperatures and has also been shown to matter when estimating long-term subsurface temperature trends.

Acknowledgements

This study was funded by the Swedish Geological Survey (project 362-1593/2013) and is part of the research within the Climate Community of the Swedish e-Science Research Centre. Partial support from Stiftelsen Lars Hiertas Minne is gratefully acknowledged. The authors thank the PFLOTRAN-ICE developers and in particular Dr Satish Karra for providing valuable technical advices, Professor Peter Jansson for providing Tarfala site data and three anonymous reviewers for comments which helped improve the study.

References

Andersland O.B. & Ladanyi B. 1994. An introduction to frozen ground engineering. Dordrecht: Springer. PubMed Abstract

Andréasson P.-G. & Gee D.G. 1989. Bedrock geology and morphology of the Tarfala area, Kebnekaise Mts., Swedish Caledonides. Geografiska Annaler Series A 71, 235–239. Publisher Full Text

Atchley A.L., Painter S.L., Harp D.R., Coon E.T., Wilson C.J., Liljedahl A.K. & Romanovsky V.E. 2015. Using field observations to inform thermal hydrology models of permafrost dynamics with ATS (v0.83). Geoscientific Model Development 8, 2701–2722. Publisher Full Text

Bear J. 1979. Hydraulics of groundwater. New York: McGraw-Hill.

Bense V.F., Ferguson G. & Kooi H. 2009. Evolution of shallow groundwater flow systems in areas of degrading permafrost. Geophysical Research Letters 36, L22401, doi: http://dx.doi.org/10.1029/2009GL039225 Publisher Full Text

Bense V.F., Kooi H., Ferguson G. & Read T. 2012. Permafrost degradation as a control on hydrogeological regime shifts in a warming climate. Journal of Geophysical Research—Earth Surface 117, F03036, doi: http://dx.doi.org/10.1029/2011JF002143 Publisher Full Text

Bojinski S., Verstraete M., Peterson T.C., Richter C., Simmons A. & Zemp M. 2014. The concept of essential climate variables in support of climate research, applications, and policy. Bulletin of the American Meteorological Society 95, 1431–1443. Publisher Full Text

Chadburn S., Burke E., Essery R., Boike J., Langer M., Heikenfeld M., Cox P. & Friedlingstein P. 2015. An improved representation of physical permafrost dynamics in the JULES land surface model. Geoscientific Model Development 8, 1493–1508. Publisher Full Text

Christiansen H.H., Etzelmüller B., Isaksen K., Juliussen H., Farbrot H., Humlum O., Johansson M., Ingeman-Nielsen T., Kristensen L., Hjort J., Holmlund P., Sannel A.B.K., Sigsgaard C., Åkerman H.J., Foged N., Blikra L.H., Pernosky M.A. & Ødegård R.S. 2010. The thermal state of permafrost in the Nordic area during the international polar year 2007–2009. Permafrost and Periglacial Processes 21, 156–181. Publisher Full Text

Coon E., Berndt M., Garimella R., Moulton J.D., Manzini G. & Painter S.L. 2013. Computational advances in the Arctic terrestrial simulator: modeling permafrost degradation in a warming Arctic. American Geophysical Union Fall Meeting 2013 Abstracts 1, p. 2195.

Dingman S.L. 2002. Physical hydrology. Upper Saddle River, NJ: Prentice-Hall. PubMed Central Full Text

Endrizzi S., Gruber S., Dall’Amico M. & Rigon R. 2014. GEOtop 2.0: simulating the combined energy and water balance at and below the land surface accounting for soil freezing, snow cover and terrain effects. Geoscientific Model Development 7, 2831–2857. Publisher Full Text

Engelhardt M., Hauck C. & Salzmann N. 2010. Influence of atmospheric forcing parameters on modelled mountain permafrost evolution. Meteorologische Zeitschrift 19, 491–500. Publisher Full Text

Etzelmüller B. 2013. Recent advances in mountain permafrost research. Permafrost and Periglacial Processes 24, 99–107. Publisher Full Text

Etzelmüller B., Farbrot H., Gudmundsson A., Humlum O., Tveito O.E. & Bjornsson H. 2007. The regional distribution of mountain permafrost in Iceland. Permafrost and Periglacial Processes 18, 185–199. Publisher Full Text

Farbrot H., Isaksen K., Etzelmüller B. & Gisnås K. 2013. Ground thermal regime and permafrost distribution under a changing climate in northern Norway. Permafrost and Periglacial Processes 24, 20–38. Publisher Full Text

Frampton A. & Destouni G. 2015. Impact of degrading permafrost on subsurface solute transport pathways and travel times. Water Resources Research 51, 7680–7701. Publisher Full Text

Frampton A., Painter S., Lyon S.W. & Destouni G. 2011. Non-isothermal, three-phase simulations of near-surface flows in a model permafrost system under seasonal variability and climate change. Journal of Hydrology 403, 352–359. Publisher Full Text

Frampton A., Painter S.L. & Destouni G. 2013. Permafrost degradation and subsurface-flow changes caused by surface warming trends. Hydrogeology Journal 21, 271–280. Publisher Full Text

Gardner W.P., Hammond G. & Lichtner P. 2015. High performance simulation of environmental tracers in heterogeneous domains. Groundwater 53, 71–80. PubMed Abstract | Publisher Full Text

Ge S., McKenzie J., Voss C. & Wu Q. 2011. Exchange of groundwater and surface-water mediated by permafrost response to seasonal and long term air temperature variation. Geophysical Research Letters 38, L14402, doi: http://dx.doi.org/10.1029/2011GL047911 Publisher Full Text

Goodrich L.E. 1982. The influence of snow cover on the ground thermal regime. Canadian Geotechnical Journal 19, 421–432. Publisher Full Text

Guymon G.L. 1976. Summer moisture–temperature for Arctic tundra. Journal of the Irrigation and Drainage Division 102, 403–411.

Hammond G.E., Lichtner P.C., Lu C. & Mills R.T. 2012. Pflotran: reactive flow & transport code for use on laptops to leadership-class supercomputers. In F. Zhang et al. (eds.): Groundwater reactive transport models. Pp. 141–159. Oak Park, IL: Bentham Science Publishers.

Hansen J., Ruedy R., Glascoe J. & Sato M. 1999. GISS analysis of surface temperature change. Journal of Geophysical Research—Atmospheres 104, 30997–31022. Publisher Full Text

Harris C., Arenson L.U., Christiansen H.H., Etzelmüller B., Frauenfelder R., Gruber S., Haeberli W., Hauck C., Hölzle M., Humlum O., Isaksen K., Kääb A., Kern-Lütschg M.A., Lehning M., Matsuoka N., Murton J.B., Nötzli J., Phillips M., Ross N., Seppälä M., Springman S.M. & Vonder Mühll D. 2009. Permafrost and climate in Europe: monitoring and modelling thermal, geomorphological and geotechnical responses. Earth-Science Reviews 92, 117–171. Publisher Full Text

Harris C., Haeberli W., Vonder Mühll D. & King L. 2001. Permafrost monitoring in the high mountains of Europe: the PACE Project in its global context. Permafrost and Periglacial Processes 12, 3–11. Publisher Full Text

Hillel D. 1980. Applications of soil physics. New York: Academic Press. PubMed Central Full Text

Hinkel K.M. & Hurd J.K. 2006. Permafrost destabilization and thermokarst following snow fence installation, Barrow, Alaska, U.S.A. Arctic, Antarctic, and Alpine Research 38, 530–539. Publisher Full Text

Hinzman L.D., Goering D.J. & Kane D.L. 1998. A distributed thermal model for calculating soil temperature profiles and depth of thaw in permafrost regions. Journal of Geophysical Research—Atmospheres 103, 28975–28991. Publisher Full Text

Isaksen K., Holmlund P., Sollid J.L. & Harris C. 2001. Three deep alpine-permafrost boreholes in Svalbard and Scandinavia. Permafrost and Periglacial Processes 12, 13–25. Publisher Full Text

Isaksen K., Mühll D.V., Gubler H., Kohl T. & Sollid J.L. 2000. Ground surface-temperature reconstruction based on data from a deep borehole in permafrost at Janssonhaugen, Svalbard. Annals of Glaciology 31, 287–294. Publisher Full Text

Isaksen K., Sollid J.L., Holmlund P. & Harris C. 2007. Recent warming of mountain permafrost in Svalbard and Scandinavia. Journal of Geophysical Research—Earth Surface 112, F02S04, doi: http://dx.doi.org/10.1029/2006JF000522 Publisher Full Text

Jonsell U., Hock R. & Duguay M. 2013. Recent air and ground temperature increases at Tarfala Research Station, Sweden. Polar Research 32, article no. 19807, doi: http://dx.doi.org/10.3402/polar.v32i0.19807 Publisher Full Text

Juliussen H. & Humlum O. 2007. Towards a TTOP ground temperature model for mountainous terrain in central–eastern Norway. Permafrost and Periglacial Processes 18, 161–184. Publisher Full Text

Kane D.L., Hinzman L.D. & Zarling J.P. 1991. Thermal response of the active layer to climatic warming in a permafrost environment. Cold Regions Science and Technology 19, 111–122. Publisher Full Text

Karra S., Painter S.L. & Lichtner P.C. 2014. Three-phase numerical model for subsurface hydrology in permafrost-affected regions. The Cryosphere Discuss 8, 149–185. Publisher Full Text

Karunaratne K.C. & Burn C.R. 2003. Freezing n-factors in discontinuous permafrost terrain, Takhini River Valley, Yukon Territory, Canada. In M. Phillips et al. (eds.): ICOP 2003. Permafrost. 8th International Conference on Permafrost. Vol. 1. Pp. 519–524. Rotterdam: Balkema.

Koven C.D., Riley W.J. & Stern A. 2013. Analysis of permafrost thermal dynamics and response to climate change in the CMIP5 earth system models. Journal of Climate 26, 1877–1900. Publisher Full Text

Kurylyk B.L., MacQuarrie K.T.B. & McKenzie J.M. 2014. Climate change impacts on groundwater and soil temperatures in cold and temperate regions: implications, mathematical theory, and emerging simulation tools. Earth-Science Reviews 138, 313–334. Publisher Full Text

Ling F. & Zhang T. 2004. A numerical model for surface energy balance and thermal regime of the active layer and permafrost containing unfrozen water. Cold Regions Science and Technology 38, 1–15. Publisher Full Text

Luetschg M., Lehning M. & Haeberli W. 2008. A sensitivity study of factors influencing warm/thin permafrost in the Swiss Alps. Journal of Glaciology 54, 696–704. Publisher Full Text

Lunardini V. 1978. Theory of n-factors and correlation of data. In: Proceedings of the Third International Conference on Permafrost, July 10–13, 1978, Edmonton, AB. Vol. 1. Pp. 40–46. Ottawa: National Research Council of Canada.

Marmy A., Salzmann N., Scherler M. & Hauck C. 2013. Permafrost model sensitivity to seasonal climatic changes and extreme events in mountainous regions. Environmental Research Letters 8, article no. 035048, doi: http://dx.doi.org/10.1088/1748-9326/8/3/035048 Publisher Full Text

McKenzie J.M., Voss C.I. & Siegel D.I. 2007. Groundwater flow with energy transport and water–ice phase change: numerical simulations, benchmarks, and application to freezing in peat bogs. Advances in Water Resources 30, 966–983. Publisher Full Text

O’Donnell J.A., Romanovsky V.E., Harden J.W. & McGuire A.D. 2009. The effect of moisture content on the thermal conductivity of moss and organic soil horizons from black spruce ecosystems in interior Alaska. Soil Science 174, 646–651. Publisher Full Text

Osterkamp T.E. 1984. Potential impact of a warmer climate on permafrost in Alaska. In J.H. McBeath (ed.): Proceedings of Conference on the Potential Effects of Carbon-induced Climatic Changes in Alaska. Miscellaneous Publication 83–1. Pp. 106–113. Fairbanks: School of Agriculture and Land Resources Management, University of Alaska. PubMed Abstract

Osterkamp T.E. & Romanovsky V.E. 1999. Evidence for warming and thawing of discontinuous permafrost in Alaska. Permafrost and Periglacial Processes 10, 17–37. Publisher Full Text

Painter S.L. 2011. Three-phase numerical model of water migration in partially frozen geological media: model formulation, validation, and applications. Computational Geosciences 15, 69–85. Publisher Full Text

Painter S.L. & Karra S. 2014. Constitutive model for unfrozen water content in subfreezing unsaturated soils. Vadose Zone Journal 13, doi: http://dx.doi.org/10.2136/vzj2013.04.0071 Publisher Full Text

Riseborough D., Shiklomanov N., Etzelmüller B., Gruber S. & Marchenko S. 2008. Recent advances in permafrost modelling. Permafrost and Periglacial Processes 19, 137–156. Publisher Full Text

Robertson E.C. 1988. Thermal properties of rocks. USGS Open File Report 88–441. Reston, VA: US Geological Survey.

Romanovsky V.E., Sazonova T.S., Balobaev V.T., Shender N.I. & Sergueev D.O. 2007. Past and recent changes in air and permafrost temperatures in eastern Siberia. Global and Planetary Change 56, 399–413. Publisher Full Text

Scherler M., Hauck C., Hoelzle M., Stähli M. & Völksch I. 2010. Meltwater infiltration into the frozen active layer at an alpine permafrost site. Permafrost and Periglacial Processes 21, 325–334. Publisher Full Text

Smith M.W. & Riseborough D.W. 1996. Permafrost monitoring and detection of climate change. Permafrost and Periglacial Processes 7, 301–309. Publisher Full Text

Sollid J.L., Holmlund P., Isaksen K. & Harris C. 2000. Deep permafrost bore-holes in western Svalbard, northern Sweden and southern Norway. Norwegian Journal of Geography 54, 186–191.

Solomon S., Qin D., Manning M., Chen Z., Marquis M., Averyt K.B., Tignor M. & Miller H.L. Jr. (eds.) 2007. Climate change 2007. The physical science basis: contribution of Working Group I to the fourth assessment report of the Intergovernmental Panel on Climate Change. Cambridge: Cambridge University Press.

Stieglitz M., Déry S.J., Romanovsky V.E. & Osterkamp T.E. 2003. The role of snow cover in the warming of Arctic permafrost. Geophysical Research Letters 30, article no. 1721, doi: http://dx.doi.org/10.1029/2003GL017337 Publisher Full Text

Turner J., Overland J.E. & Walsh J.E. 2007. An Arctic and Antarctic perspective on recent climate change. International Journal of Climatology 27, 277–293. Publisher Full Text

Tutolo B.M., Kong X.-Z., Seyfried W.E. Jr. & Saar M.O. 2015. High performance reactive transport simulations examining the effects of thermal, hydraulic, and chemical (THC) gradients on fluid injectivity at carbonate CCUS reservoir scales. International Journal of Greenhouse Gas Control 39, 285–301. Publisher Full Text

van Genuchten M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44, 892–898. Publisher Full Text

Vonder Mühll D. & Haeberli W. 1990. Thermal characteristics of the permafrost within an active rock glacier (Murtel/Corvatsch, Grisons, Swiss Alps). Journal of Glaciology 36, 151–158.

Williams P.J. & Smith M.W. 1991. The frozen earth. Fundamentals of geocryology. Cambridge: Cambridge University Press.

Woo M.-K. 2012. Permafrost hydrology. Dordrecht: Springer. Publisher Full Text

Zhang T. 2005. Influence of the seasonal snow cover on the ground thermal regime: an overview. Reviews of Geophysics 43, RG4002, doi: http://dx.doi.org/10.1029/2004RG000157 Publisher Full Text

Zhang T., Barry R.G., Knowles K., Ling F. & Armstrong R.L. 2003. Distribution of seasonally and perennially frozen ground in the Northern Hemisphere. In M. Phillips et al. (eds.): ICOP 2003. Permafrost. 8th International Conference on Permafrost. Vol. 2. Pp. 1289–1294. Rotterdam: Balkema.

Zhang T., Osterkamp T.E. & Stamnes K. 1996. Influence of the depth hoar layer of the seasonal snow cover on the ground thermal regime. Water Resources Research 32, 2075–2086. Publisher Full Text

Zhang T., Osterkamp T.E. & Stamnes K. 1997. Effects of climate on the active layer and permafrost on the North Slope of Alaska, USA. Permafrost and Periglacial Processes 8, 45–67. Publisher Full Text