# Evaluating the performance of coupled snow-soil models in SURFEXv8 to simulate the permafrost thermal regime at a high Arctic site

Mathieu Barrere, Florent Domine, Bertrand Decharme, Samuel Morin, Vincent Vionnet, Matthieu Lafaysse

► To cite this version:

Mathieu Barrere, Florent Domine, Bertrand Decharme, Samuel Morin, Vincent Vionnet, et al.. Evaluating the performance of coupled snow-soil models in SURFEXv8 to simulate the permafrost thermal regime at a high Arctic site. Geoscientific Model Development Discussions, 2017, 10 (9), pp.3461-3479. 10.5194/gmd-10-3461-2017 . insu-03326445

HAL Id: insu-03326445

https://insu.hal.science/insu-03326445

Submitted on 26 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Evaluating the performance of coupled snow-soil models in SURFEXv8 to simulate the permafrost thermal regime at a high Arctic site

Mathieu Barrere<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" />,<fnref n="5" />, Florent Domine<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="6" />, Bertrand Decharme<fnref n="7" />, Samuel Morin<fnref n="4" />, Vincent Vionnet<fnref n="4" />, and Matthieu Lafaysse<fnref n="4" />

*Correspondence to:* Florent Domine (florent.domine@gmail.com)

Received: 1 March 2017 – Discussion started: 21 March 2017

Revised: 11 August 2017 – Accepted: 18 August 2017 – Published: 21 September 2017

**Abstract.** Climate change projections still suffer from a limited representation of the permafrost-carbon feedback. Predicting the response of permafrost temperature to climate change requires accurate simulations of Arctic snow and soil properties. This study assesses the capacity of the coupled land surface and snow models ISBA-Crocus and ISBA-ES to simulate snow and soil properties at Bylot Island, a high Arctic site. Field measurements complemented with ERA-Interim reanalyses were used to drive the models and to evaluate simulation outputs. Snow height, density, temperature, thermal conductivity and thermal insulation are examined to determine the critical variables involved in the soil and snow thermal regime. Simulated soil properties are compared to measurements of thermal conductivity, temperature and water content. The simulated snow density profiles are unrealistic, which is most likely caused by the lack of representation in snow models of the upward water vapor fluxes generated by the strong temperature gradients within the snowpack. The resulting vertical profiles of thermal conductivity are inverted compared to observations, with high simulated values at the bottom of the snowpack. Still, ISBA-Crocus manages to successfully simulate the soil temperature in winter. Results are satisfactory in summer, but the temperature of the top soil could be better reproduced by adequately represent-

ing surface organic layers, i.e., mosses and litter, and in particular their water retention capacity. Transition periods (soil freezing and thawing) are the least well reproduced because the high basal snow thermal conductivity induces an excessively rapid heat transfer between the soil and the snow in simulations. Hence, global climate models should carefully consider Arctic snow thermal properties, and especially the thermal conductivity of the basal snow layer, to perform accurate predictions of the permafrost evolution under climate change.

## 1 Introduction

The Arctic is warming at twice the average planetary rate (Sweet et al., 2015). As climatological, hydrological and biological systems are fully coupled (Hinzman et al., 2005), climate change affects each Arctic ecosystem (Post et al., 2009; Serreze et al., 2000). Consequences are already observed, such as reduced sea ice extent and snow cover duration (Serreze et al., 2000), shifts in vegetation (Pearson et al., 2013), permafrost degradation (Smith et al., 2010) and faunal redistribution (Post and Forchhammer, 2008).

<footnote n="1">¹Centre d'Études Nordiques, Université Laval, Québec City, QC, Canada</footnote>

<footnote n="2">²Department of Geography, Université Laval, Québec City, QC, Canada</footnote>

<footnote n="3">³Takuvik Joint International Laboratory, Université Laval (Canada) and CNRS-INSU (France), Québec City, QC, Canada</footnote>

<footnote n="4">⁴Météo-France – CNRS, CNRM UMR 3589, CEN, Grenoble, France</footnote>

<footnote n="5">⁵Univ. Grenoble Alpes – CNRS – IRD, IGE, Grenoble, France</footnote>

<footnote n="6">⁶Department of Chemistry, Université Laval, Québec City, QC, Canada</footnote>

<footnote n="7">⁷Météo-France – CNRS, CNRM UMR 3589, Toulouse, France</footnote>
Permafrost degradation is of major concern because of its feedback on the global climate system. Indeed, large amounts of organic carbon are stored in perennially frozen soils because of the limited microbial activity (Jonasson et al., 2001). Recent studies estimate that about 1300 Pg of soil organic carbon (SOC) are stored in permafrost (Hugelius et al., 2014), constituting one of the largest terrestrial carbon pools. As permafrost warms and thawing occurs, SOC becomes available for microbial mineralization, resulting in the release of potentially very significant amounts of greenhouse gases to the atmosphere (Elberling et al., 2013; Schuur et al., 2015). This effect is considered as one of the strongest positive climate feedbacks and needs to be taken into account in global temperature predictions. However, the permafrost-carbon feedback has not been included in the climate projections of the IPCC Fifth Assessment Report (Schaefer et al., 2014), so that current warming projections may be significantly underestimated.

Considerable uncertainties remain in the SOC decomposition rate and associated greenhouse gas emissions in this climate change context. One of the main reasons is that the rate and extent of permafrost thaw is not well quantified, preventing accurate estimates of the potential magnitude of the permafrost-carbon feedback. In particular, how the permafrost thermal regime will respond to climate change is still poorly represented because of its high sensitivity to the properties of the surface. The snow cover and the vegetation type, by modifying surface energy exchanges, are the main local factors influencing the permafrost thermal regime (Sturm et al., 2001a; Zhang, 2005). A snow cover acts as a thermal insulator by limiting soil winter cooling, but its insulating properties are highly variable and insufficiently detailed in global climate models. These characteristics, and how they will change with climate, therefore need to be considered to successfully simulate the current and future permafrost thermal regime.

We attempt to contribute to these aspects by investigating the capacity of the soil multilayer version of the land surface model ISBA (Interactions Soil Biosphere Atmosphere; Decharme et al., 2013) coupled with the detailed snowpack scheme Crocus (Vionnet et al., 2012) or to the simpler scheme ES (Explicit Snow; Boone and Etchevers, 2001; Decharme et al., 2016) to simulate snow and soil properties at a high Arctic location. Crocus is currently the most detailed snowpack scheme coupled to ISBA, and corresponds to the highest existing level of complexity for snowpack models. However, its relatively large computation time limits its use in global climate models. ES is a multilayer snowpack scheme of intermediate complexity. Based on Crocus parameterizations, ES will be used for the upcoming CMIP6 CNRM-CM simulations (Eyring et al., 2016). Similar studies were performed over Siberian regions, using ISBA-Crocus (Brun et al., 2013) and ISBA-ES (Decharme et al., 2016). Models were evaluated based on snow height, snow water

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Location of the study site in the southwest plain of Bylot Island in the Canadian Arctic archipelago.</figcaption></figure>

equivalent and soil temperature, but processes related to the simulated snowpacks have not been studied in detail yet.

Here, simulations are tested using field data obtained at one specific site on Bylot Island (73°N, 80°W), located north of Baffin Island in the Canadian Arctic archipelago. The mean annual air temperature for the period 1998–2013 is –14.3 °C and permafrost thickness has been estimated to be over 400 m (Fortier and Allard, 2004). The site description and instrumentation are given in Sect. 2, after which the main features of the models are detailed. Simulation results are shown and compared to field observations in Sect. 3. A sensitivity analysis is performed to determine processes critical to simulate Arctic snow and soil, which are discussed in Sect. 4 of the paper.

# 2 Methods

## 2.1 Site description

The study area is in the Qarlikturvik valley of the southwest plain of Bylot Island, around 73°10′N, 80°00′W (Fig. 1). The glacial retreat around 6000 years ago left fine-grained wind deposits and organic sediments to form the soil (Allard, 1996). The valley bottom consists of wetlands with typical permafrost landforms, including tundra polygons, thaw lakes and ponds. Vegetation is mainly comprised of sedges, graminoids and mosses (Cadieux and Gauthier, 2008).

Our actual study site is located in wetlands of the valley floor, at 73°09′01.4″N, 80°00′16.6″W, in the middle of a low-center tundra polygon (Domine et al., 2016a). Vegetation consists of a typical herb tundra environment, covered by sedges, graminoids (*Poa* sp.), mosses (*Aulacumnium turgidum*, 2–3 cm thick) and small prostrate Arctic willows (*Salix herbacea* and *S. arctica*). Albedo measurements were performed in July 2015, using a SVC HR-1024 spectroradiometer.
diometer. Averaged over the 346–2513 nm spectral range, the surface albedo of our very site was 0.17 (M. Belke-Brea, personal communication, 2015). The soil granulometry was analyzed using a laser scattering particle size analyzer, which resulted in a combination of fine grain deposits comprised of 64 % silt and 36 % fine sand (Domine et al., 2016a), based on the United States Department of Agriculture classification. An organic litter layer (3.5 to 6 cm thick) was present at the interface between the surface vegetation and the ground. Because of the difficulty to determine the boundaries between the living moss, the litter and the underlying mineral soil, we estimated a total thickness of about 10 cm for both moss and litter. Field observations and simulations presented in this paper focus on this very spot.

## 2.2 Site instrumentation and data

As we investigate the coupled evolution of the atmosphere, snow and soil, various instruments have been installed to monitor meteorological conditions, along with snowpack and soil properties. The set of instruments was installed on the same tundra polygon for the comprehensive monitoring of the spot. The protocol of instrumentation and measurement has been fully discussed in Domine et al. (2016a, b), thus only a brief description is given here.

Several automatic weather stations are operating on Bylot Island (CEN, 2016). In particular, a 10 m tower was installed in summer 2004 close to our study site, at 73°09′07.9″ N, 79°59′19.0″ W, to monitor air temperature and relative humidity, wind speed and direction, and snow height. During summer 2013, an automatic weather station (hereafter referred to as BylSta) was installed on the polygon of our study site (Domine et al., 2016a). It measures air temperature and relative humidity with a ventilated HC2S3 sensor from Rotronic, wind speed with a cup anemometer, surface temperature from an infrared IR120 sensor, and upwelling and downwelling shortwave and longwave radiation with a CNR4 radiometer associated with a CNF4 heating/ventilating unit from Kipp & Zonen. An issue with the HC2S3 sensor caused erroneous air temperature measurements in the first year. This was fixed in summer 2014, providing a fairly complete meteorological dataset since then.

Snow height was automatically monitored with a SR50A acoustic gauge installed on the automatic weather station. A few meters further, three TP08 heated needle probes (NPs) from Hukseflux were placed at 7, 17 and 27 cm above the ground in summer 2013 to measure the snow temperature and its thermal conductivity $k_{\text{snow}}$. They were lowered to 2, 12 and 22 cm in July 2014 to better match the snowpack structure observed in May 2014. Operating methods and data analysis pertaining to the NPs are detailed in Domine et al. (2015, 2016a) and Morin et al. (2010). Applied to snow, the NP method is suspected of presenting a low systematic error of 20 % on average, related to the granular structure of the medium (Calonne et al., 2011; Domine et al., 2015;

Riche and Schneebeli, 2013). The anisotropy of the snow structure is another possible source of error, because horizontally inserted NPs measure a mixture of vertical and horizontal thermal conductivities. As heat exchanges between the ground and the atmosphere occur in the vertical dimension through the snow, the anisotropy of the snow thermal conductivity can produce errors up to 20 % in measurements, resulting in a maximum total error of 29 % (Domine et al., 2015). Compared to the large range of $k_{\text{snow}}$ values (0.025–0.7 W m⁻¹ K⁻¹), and given that this method is the only suitable solution for remote field work, errors related to the use of NP are acceptable.

Field campaigns took place in May 2014 and 2015 at the end of the snow season. The snow accumulation was highly variable because of wind effects and microtopography (Liston and Sturm, 2002). The SR50A automatic point measurements are therefore not necessarily representative of the average snow conditions. To explore the spatial variability of snow properties, we performed hundreds of snow height measurements covering different areas with an avalanche probe. In addition, snow pits were dug at a dozen specific sites to describe the stratigraphy and to measure vertical profiles of density with a 100 cm³ box cutter, and temperature and thermal conductivity with a TP02 heated needle. As we are focusing on the station site, we will only use data from the three snow pits dug within 20 m of the station both years.

Soil temperature and volumetric water content (VWC) have been monitored since July 2013 with Decagon 5TM probes, installed at depths of 2, 5, 10 and 15 cm. They were not calibrated for our specific soil, so we used the manufacturer's calibration for mineral soils which may produce an error in water content of 3 %. The temperature sensor accuracy is within ±1 °C. The active layer was 17 cm thick at the time of installation so the sensors could not be placed deeper. A TP08 heated needle probe was inserted at 10 cm depth, just below the litter layer, to automatically monitor the soil thermal conductivity $k_{\text{soil}}$. The method used is the same as for the $k_{\text{snow}}$ measurements, and data analysis is also detailed in Domine et al. (2016a). With our instrumental methods, we must use the same heating power in the NPs for snow and soil. The heating power was optimized for snow to minimize heating and hence perturbation to metamorphism. Since soils have a higher thermal conductivity than snow, especially when frozen, the thermal signal of the NP is low and noise for frozen soil thermal conductivity data is higher than what could be achieved using a higher heating power. Two field campaigns were conducted in July 2014 and 2015 during which we measured soil water content profiles using an EC-5 sensor from Decagon, and temperature and thermal conductivity with a TP02 heated needle at several dozen sites (Domine et al., 2016a).
## 2.3 Simulations

Simulations were performed using SURFEX v8 (SURFace EXternalisée), the surface modelling platform developed by Météo-France and partners (Masson et al., 2013). Used in stand-alone mode, it takes meteorological data as driving input. The snow cover and the underlying soil are coupled through a semi-implicit scheme (Decharme et al., 2016; Vionnet et al., 2012). Different configurations are tested and compared with field observations to identify critical processes affecting the permafrost thermal regime. Based on a series of incremental runs, we particularly focus on the following model features: surface litter, soil organic carbon and density of the drifting snow.

### 2.3.1 Meteorological driving data

To calculate the energy and mass budget of the surface, the model needs the following input data: air temperature and specific humidity, wind speed, incoming shortwave and long-wave radiation, precipitation rate (solid and liquid) and atmospheric pressure. We used observed local meteorological data when available, and missing data are filled with ERA-Interim (ERAi) reanalysis (Dee et al., 2011). Available from 1979 to present, for a 0.7° grid with a 3-hourly time resolution, ERAi provides a continuous meteorological dataset available globally.

Air temperature and wind speed observations are available since 2004 from the 10 m tower. After July 2014, we used the ventilated air temperature and humidity, and the wind speed measured at BylSta. Radiation measurements were not used because the radiometer shifted by about 5° from its horizontal position when the tripod that supported it sank with the active layer thawing, causing errors in the data. The atmospheric pressure and precipitation were not measured, so we used ERAi reanalysis data for missing variables and possible data gaps.

To make the ERAi data consistent with the original field data, they were corrected following the method of Vuichard and Papale (2015). This method consists of calculating the linear regression between ERAi data and available field measurements. The regression coefficients (slope and intercept) are used to correct systematic biases in ERAi data, in order to better match the local meteorological conditions. Hence, we found a mean bias of -2.3 °C in air temperature, -2.7 × 10⁻⁴ g kg⁻¹ in air humidity and +0.7 m s⁻¹ in wind speed given by ERAi. The correction helped to reduce the standard deviation between the corrected ERAi and the observational data by respectively 20, 4.6 and 10%. Differences in the correction performance mainly reflect the internal variability of the in situ data.

The precipitation phase was recalculated using local and ERAi-corrected temperatures, with a threshold set at 1 °C. During the ERAi reanalysis process, the precipitation fluxes are only predicted by the meteorological model with no as-

simulation of observations (Brun et al., 2013). Furthermore, our study site is surrounded by hills, altering the local aerology and therefore the precipitation amount. With no reliable precipitation gauge to evaluate the local variability, we consider the precipitation data as the highest uncertainty source of all the forcing data. Preliminary results indicated that an accurate snow height was critical to simulate correctly the soil thermal regime. Hence, we arbitrarily changed the ERAi precipitation data for Crocus to match the observed snow height at snow pits dug in the immediate vicinity of BylSta, which was achieved after reducing the solid precipitation rate by 30% in winter 2013–2014. Neither liquid precipitations nor snowfall were modified for the winter 2014–2015. ERAi radiation data were kept unchanged for lack of reliable measured values, as well as atmospheric pressure values which we assume suffer little from local variability.

### 2.3.2 Soil scheme

The land-surface parameterization is managed by the ISBA scheme (Noilhan and Mahfouf, 1996; Noilhan and Planton, 1989). It describes the exchanges of energy and water between the atmosphere, the vegetation and the soil, by solving the 1-D Fourier's law for heat transfer and the mixed-form of Richards' equation for water mass transfer within the soil (Decharme et al., 2011). They are computed using 20 layers down to a depth of 12 m. The main parameters are the soil texture and the vegetation type, from which other parameters can be derived, e.g., the soil porosity, the saturated matric potential and the saturated hydraulic conductivity. Following our measurements, the soil texture is set to 36% sand and 64% silt. Boreal grass covers the surface with a root depth of 20 cm, and the snow-free albedo is 0.17.

The soil freezing-thawing processes are critical in permafrost regions. They are handled using the drying-wetting analogy based on the Gibbs free-energy method (Boone et al., 2000). It allows the calculation of the temperature for phase changes as a function of the soil matric potential, which depends on porosity and water content. During phase changes, the liquid water content decreases (increases) correspondingly to the increase (decrease) in ice content, thereby conserving the total water content in each soil layer. ISBA calculates the hydrology of the entire soil column in order to accurately represent the permafrost characteristics. Therefore, the model needs long simulation periods to guarantee that the water and heat profiles are equilibrated over the 12 m soil depth. Water infiltration is limited by the presence of ice-rich permafrost and because of the topographic hollow of our study site water is at times poorly drained and can occasionally fill the center of the polygon. This process is reproduced by disabling lateral runoff when surface water is in excess (run base and following).

The soil thermal properties, i.e., thermal conductivity and heat capacity, are computed as a combination of water, ice and soil properties, volumetric water content and soil
porosity, following the parameterizations of Peters-Lidard et al. (1998). Hence, the soil thermal conductivity is expressed as a function of its saturation, porosity, quartz content, dry soil conductivity and phase of water (frozen or unfrozen), where the ice, water and quartz thermal conductivities are respectively 2.2, 0.57 and 7.7 W m⁻¹ K⁻¹. Surface organic material like moss or litter are known to greatly affect the thermal and hydraulic properties of the soil (Hinzman et al., 1991). ISBA cannot handle a surface moss cover, but it has the capacity to simulate a litter by conferring to the uppermost soil layers the thermal properties of organic matter. The thermal conductivity of organic matter is set to 0.05 W m⁻¹ K⁻¹ when dry, and to 0.25 W m⁻¹ K⁻¹ when wet. We evaluated the model sensitivity to the inclusion of a 10 cm thick litter which affects thermal properties of the first 10 cm below the surface (run litter and following). Based on our measurements, the litter appears to have a high insulating capacity in summer, while it becomes negligible in winter. To better reproduce the observations, the litter effect is therefore disabled as soon as 2 cm of snow covers the surface. The soil hydrology is not affected in the model, while in reality organic layers have a high hydraulic conductivity and a greater infiltration rate than mineral soils (Hinzman et al., 1991). The impact of this simplification will be evaluated with the results.

ISBA includes the dependency of soil hydraulic and thermal properties on SOC content as fully described in Decharme et al. (2016). Briefly, depending on its content, which decreases sharply with soil depth, including organic carbon reduces the dry soil thermal conductivity, increases its porosity and therefore increases its saturated hydraulic conductivity. Accounting for the effect of SOC could significantly help improve the simulation of the soil thermal regime, especially in Arctic areas where soils store large amounts of SOC (Hugelius et al., 2014). From the Harmonized World Soil Database (HWSD; FAO, 2012) with a 1 km resolution, the SOC content for our study site is estimated at 5.05 kg m⁻³ for the first soil horizon (0–30 cm), and 4.34 kg m⁻³ from 30 cm to 1 m depth. A few soil samples were collected in July 2013 on the southern plain of Bylot Island for a radio-carbon analysis (ADAPT 2014). From two humid and three mesic sites, the carbon concentration was 5.17 kg m⁻³ averaged over the first 30 cm of the soils. Given the excellent agreement with ADAPT data, we have a good confidence in the HWSD estimations used in our simulations (run SOC and following).

### 2.3.3 Snowpack model: Crocus

Crocus is a multilayer physical snow scheme designed to simulate the evolution of the snow cover as a function of energy and mass transfer between the snowpack, the atmosphere and the ground (Brun et al., 1989; Vionnet et al., 2012). Numerical snow layers are handled dynamically by the model, in order to keep their total number below a given

number (typically 50), while respecting as much as possible the internal structure of the snowpack. Numerical snow layers are characterized by their thickness, density, temperature, liquid water content and four variables representing the microstructural properties (Vionnet et al., 2012). The freshly fallen snow is usually considered as dendritic, and evolves toward non-dendritic snow under the action of internal processes such as snow metamorphism, compaction, thermal diffusion and phase change.

These processes also affect the density of the snow layers, which is a key variable controlling other snow physical properties in Crocus. Density increases with compaction caused by the weight of upper layers, depending on the viscosity of the layer. The compaction rate is faster in the presence of liquid water, and it also depends on snow grain type. Angular grains such as depth hoar compact less. Wind events can also compact the surface layers, and due to frequent wind in the Arctic large amounts of snow are transported, compacted and sublimated (Liston and Sturm 2002; Sturm et al., 2001a, b). Blowing snow occurs in the simulation when the wind speed is above a threshold value which depends on surface snow properties: microstructure and density (Guyomarc'h and Merindol, 1998; Vionnet et al., 2012). On average, blowing snow is observed for wind speeds greater than 5 m s⁻¹ (Sturm et al., 2001b; Vionnet et al., 2013). From the 10 m tower on Bylot Island, wind speeds greater than 5 m s⁻¹ occurred during 7% of the 2013–2014 snow season, and during 6% of the 2014–2015 one. Crocus takes into account the compaction and the microstructural evolution of the surface snow caused by blowing snow events, depending on wind speed, snow density and its microstructural properties (Vionnet et al., 2012). In addition, Crocus calculates the snow mass lost by sublimation during blowing snow events (Brun et al., 2013; Gordon et al., 2006). However, Crocus does not handle the snow redistribution since the model is one-dimensional (Brun et al., 2013). In the standard version of Crocus, wind-compact ed snow layers can reach a maximum density of 350 kg m⁻³. As we observed snow densities up to 450 kg m⁻³ on Bylot Island, and values up to 600 kg m⁻³ can be found in the literature (e.g., Sturm et al., 1997; Zhang, 2005), we increased this value to 600 kg m⁻³ (run wind and following). Because the mobility of snow layers depends on their density, increasing the density of the snow can limit its driftability.

The snow thermal conductivity $k_{\text{snow}}$, in W m⁻¹ K⁻¹, which is used to solve the thermal diffusion equation in the snowpack, is calculated from the density using the equation of Yen (1981):

$$k_{\text{snow}} = k_{\text{ice}} \left( \frac{\rho}{\rho_{\text{w}}} \right)^{1.88}, \qquad (1)$$

with $k_{\text{ice}}$ the thermal conductivity of ice (2.22 W m⁻¹ K⁻¹), $\rho$ the density of snow and $\rho_{\text{w}}$ the density of liquid water (1000 kg m⁻³).
<table><caption>**Table 1.** Configurations of the numerical experiments, all realized with Crocus except for ES. Changes are incremented in addition to the previous configuration.</caption><thead><tr><td>Run name</td><td>Run configuration</td></tr></thead><tbody><tr><td>base</td><td>base simulation using SURFEX version 8.0 (rev. 4006), for a mineral soil</td></tr><tr><td>litter</td><td>addition of a 10 cm thick surface litter</td></tr><tr><td>SOC</td><td>addition of organic carbon within the soil profile</td></tr><tr><td>wind</td><td>increase the maximum density for wind-induced snow compaction, from 350 to 600 kg m<sup>-3</sup></td></tr><tr><td>ES</td><td>use the ES snow scheme instead of Crocus, same configuration as wind</td></tr></tbody></table>

Snow albedo depends on the snow microstructure, amount of light-absorbing impurities and on the solar zenith angle (SZA; Warren, 1982). Crocus computes the snow albedo by considering microstructure to account for physical aspects and the age of the surface snow layer to account for chemical aspects, as exposed snow is subjected to impurity deposition. The incoming radiation is then transmitted and absorbed within the snowpack, following the exponential decay of radiation with depth as a function of the grain size and snow density (Vionnet et al., 2012). Because incoming radiation is preferentially scattered forward in the snow, light coming from a low sun angle penetrates less deep in the snowpack and the resulting albedo is higher. Despite its importance in Arctic regions, this last effect is not simulated in Crocus because it does not account for the SZA in the albedo calculation. The physically based radiative transfer model TARTES has been implemented in Crocus, making use of the SZA in the albedo calculation (Charrois et al., 2016; Libois et al., 2015) but also requiring knowledge of the vertical profile of light absorbing impurities (Tuzet et al., 2017). However, the lack of data on snow impurities (nature, deposition, light-absorbing spectroscopy) prevented us from using TARTES on Bylot Island simulations.

### 2.3.4 Snowpack model: Explicit Snow

Explicit Snow (ES) is a multilayer snowpack scheme of intermediate complexity (Boone and Etchevers, 2001; Decharme et al., 2016). It is based on parameterizations used in Crocus, which are simplified to reduce the computation time and to facilitate its integration in global climate models. The main differences with Crocus are a constant number of layers (usually 12), blowing snow is not sublimated and the snow microstructural properties are not simulated. Hence, parameterizations of snow albedo and compaction rate are functions of the layers density only. However, as in Crocus, the evolution of snow density in each layer is due to snow compaction resulting from changes in snow viscosity and wind-induced densification of near-surface snow layers (Decharme et al., 2016).

The $k_{\text{snow}}$ calculation method is different from Crocus. In ES, $k_{\text{snow}}$ is computed from the density, based on the Yen's equation (Eq. 1), with an additional term to account for latent heat transfer through sublimation–condensation pro-

cesses during metamorphism (Sun et al., 1999):

$$k_{\text{snow}} = k_{\text{ice}} \left(\frac{\rho}{\rho_w}\right)^{1.88} + \frac{P_0}{P_a} \times \max\left(0, k_1 - \frac{k_2}{T_s - k_3}\right), \quad (2)$$

where $P_a$ (Pa) is the air pressure, $P_0$ a reference pressure equal to 1000 hPa, $T_s$ (K) the temperature of the snow layer, and the coefficients $k_1 = -0.06023 \text{ W m}^{-1} \text{ K}^{-1}$, $k_2 = 2.5425 \text{ W m}^{-1}$ and $k_3 = 289.99 \text{ K}$.

Even though ES is not based on the explicit representation of the snow microstructure, it remains one of the most sophisticated snowpack models that will be used for the upcoming 6th edition of the Coupled Model Intercomparison Project (CMIP6; Eyring et al., 2016). The results of our simulations with ISBA-ES (run ES) are expected to reflect the capacity of the latest generation of global climate models to simulate Arctic snow and soil properties.

### 2.3.5 Numerical experiments

A first simulation was run from August 1979 to August 2012, constituting 33 years of initialization. This operation allowed reaching equilibrium between the ground properties and the local climate conditions, with the base configuration. Then, using the equilibrated soil profile as the initial state, we performed a series of runs with incremental complexities from August 2012 to June 2015 (Table 1). The run wind integrates all our changes, it consists of the most detailed configuration tested here with Crocus. ES uses the same configuration as the run wind.

## 2.4 Evaluation metrics

Field observations of snow and soil properties are available from August 2013 to June 2015, allowing the evaluation of model performance to simulate two winters and one entire summer.

Simulated snow properties are compared to measurements of snow height, density, thermal conductivity and temperature. As we assess the ability of the model to reproduce the soil thermal regime, we particularly focus on the snow thermal properties. Snow is thermally characterized by its thermal conductivity $k_{\text{snow}}$. We also rely on an alternate variable that characterizes the thermal properties of the whole snowpack rather than those of each layer. The thermal insulation
<figure><img src="image_4.png" /><figcaption>Figure 2. Snow height evolution over the 2 years of observations. Results of several hundred measurements using an avalanche probe (random measurements) and of a few snow pit data points close to our site, performed in May 2014 and 2015.</figcaption></figure>

of the snowpack $R_T$ (in m$^2$ K W$^{-1}$) depends on the thickness $h$ and the thermal conductivity of each layer $i$:

$$R_T = \sum_i \frac{h_i}{k_{\text{snow},i}}. \qquad (3)$$

The simulated $R_T$ is computed from the respective density-$k_{\text{snow}}$ correlation used in the snowpack model (Eqs. 1 and 2) and simulated layer thickness, while observed $R_T$ is directly based on measurements of $k_{\text{snow}}$ and layer thickness. To determine the respective contributions of the snowpack thickness and the thermal conductivity on $R_T$, it is also necessary to look at the resulting mean thermal conductivity of the snowpack calculated as the following:

$$\overline{k_{\text{snow}}} = \frac{h}{R_T}. \qquad (4)$$

ISBA simulation results are compared to measurements of soil temperature as well as the thermal conductivity and water content, which govern the soil thermal regime. Vertical profiles show the vertical soil stratification. Performances of each run are evaluated by comparing their deviations with time series observations at different depths, using the squared correlation coefficient ($r^2$), bias (Eq. 5) and RMSE (Eq. 6) statistical errors.

$$\text{bias} = \frac{1}{n} \sum_{i}^{n} (\text{sim}_i - \text{obs}_i), \qquad (5)$$

$$\text{RMSE} = \sqrt{\frac{1}{n} \sum_{i}^{n} (\text{sim}_i - \text{obs}_i)^2}, \qquad (6)$$

where sim and obs refer to simulated and observed values respectively.

# 3 Results

## 3.1 Simulations of snow properties and comparison with field data

### 3.1.1 Snow height

Figure 2 shows snow height automatically measured with the SR50A snow gauge, along with manual observations performed in May 2014 and 2015 and simulation results. SR50A data are missing between 14 May and 26 July 2014 because of an issue with the sensor. Random snow height measurements were performed within 200 m of our site on 14 and 16 May 2014, and on 12 May 2015. They evidenced the large spatial variability in snow height, as already indicated in Table 2 of Domine et al. (2016b). Snow-free areas were frequently encountered on the polygons edges, whereas snow heights up to 60 cm were found in the center of the polygons. Thus, one-point snow height monitoring may not be representative of the average snow accumulation. On 14 May 2014, the SR50A indicated 13 cm of snow, while we obtained a mean of 16.2 cm with a standard deviation of 13.7 cm from more than 300 measurements over the whole polygons area. We observed a noticeably lower snow accumulation under the station than a few meters away, which is confirmed by the 21 cm of snow averaged on the three snow pits dug in the vicinity of the station. On 12 May 2015, while the snow gauge indicated 35 cm, random measurements indicated 25.3 ± 13.1 cm of snow. From the three snow pits made, we obtained a mean snow height of 33 cm. Based on our observations, we explain these differences by a low snow year during 2013–2014 and a high wind redistribution (Domine et al., 2016b). Since a central objective of our simulations is to reproduce the ground thermal regime as mean
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Stratigraphies and vertical profiles of density measured on 14 May 2014 (a) and 12 May 2015 (b), and simulated densities on 14 May 2014 and 6 May 2015.</figcaption></figure>

sured by our instruments, and because the large spatial vari-
ability in snow height highly affects the ground temperature
at the meter scale (Gisnas et al., 2016), the relevant snow
height at our specific site is most likely given by the snow
pits mean.

After reducing the ERAi snowfall amount by 30% in winter 2013–2014, the resulting simulated snow height was in good agreement with the snow pit means obtained in May 2014 and 2015, respectively 23.0 and 30.4 cm from the run wind and 22.4 and 32.8 cm with ES. Winter variations are well reproduced, but the beginning of the snow season occurs too early in 2013: on 18 September in simulations, while it was detected on 12 October from measurements. At this date the simulated snowpack was already 4 cm thick, which partly explains the slight overestimation compared to the mean height of May 2014 snow pits. ES simulates lower snow heights than Crocus, which could be caused by a faster compaction rate in ES compared to Crocus. In May 2015, snowmelt occurs earlier in simulations than deduced from the snow gauge. A detailed analysis of the meteorological data and simulation outputs reveals positive air temperature between 6 and 11 May 2015, triggering the partial melting of the upper snow layers in simulations. With two strong wind events (peaks exceeding 10 m s⁻¹), which occurred on 5 and 12 May, these conditions caused the decrease in the simulated snow height along with the increase in the upper layers density. The air temperature cooled back down after this event, until 18 May when the final snowmelt started. We did not observe any signs of spring melt before 18 May during the field campaign. To test whether the inexact melt onset date simulated was caused by the lack of SZA consideration, we briefly tested Crocus coupled to TARTES, which includes treatment of SZA. With TARTES, the melt onset date was accurately simulated, which leads us to suggest that

not accounting for SZA is the main cause of the inadequate
melt onset date simulation. However, a full implementation
of TARTES in Crocus for this study would have required
data on snow impurities. Hence, for further comparisons with
snow properties measured in May 2015, we will use simula-
tion results of 6 May, just before the partial melting.

**3.1.2 Snow density and thermal conductivity**

Figure 3 shows examples of snow stratigraphies observed close to the monitoring site in May 2014 and 2015. They are typical of Arctic snowpacks (Domine et al., 2016b), comprised mainly of a basal depth hoar layer (about 5 cm thick) and a top wind slab. The particularity in 2015 is the depth hoar layer which was indurated at the bottom of the snowpack, resulting from the transformation of a melt-freeze layer into depth hoar under a very high-temperature gradient (Domine et al., 2016a, b). Indurated depth hoar is harder and retains a higher cohesion than typical depth hoar, but its development usually goes along with a decrease in density and thermal conductivity (Domine et al., 2012).

Associated measured vertical profiles of density are also
shown. As expected from typical Arctic snowpacks, the
stratigraphies observed exhibit low density values for the bot-
tom depth hoar layer (between 150 and 200 kg m⁻³), and
high values for the upper wind slabs (exceeding 400 kg m⁻³).
Crocus and ES simulate inverted profiles, with generally de-
creasing densities from the bottom to the top of the snow-
pack. This is particularly visible in 2015, with the highest
density simulated at the bottom of the snowpack. In 2014
the density of the basal layer (below 5 cm) is successfully re-
produced by Crocus, but it is overestimated by ES. In 2015,
the density of the indurated depth hoar basal layer is greatly
overestimated by both Crocus and ES, respectively 1.4 and
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Vertical profiles of snow thermal conductivities measured on 14 May 2014 (a) and 12 May 2015 (b), and simulated on 14 May 2014 and 6 May 2015.</figcaption></figure>

2 times higher. It confirms that the snow compaction rate in ES is inappropriate to simulate a bottom depth hoar layer of low density. However, this partially compensates for the underestimation of density in the upper layers. Wind slab densities in upper layers are highly underestimated in simulations (between 1.5 and 2 times lower), even for the run *wind*, which allows higher maximum densities for drifting snow. In May 2014, the *wind* experiment simulates yet lower densities in sub-surface snow layers than the other Crocus runs. This is attributed to early melting, which occurs slower in this run because of colder snow temperatures compared to previous experiments. Otherwise, increasing the maximum density reached by drifting snow helps to reduce the underestimation in upper layers. Still, mean snowpack densities are lower in simulations than in observations. ES simulates higher mean densities than Crocus because the snowpack is overall more compacted, due to the faster compaction rate. Therefore, the density of the basal layer is greatly overestimated, but the density of the upper layers are less underestimated.

Since thermal conductivity is totally (in Crocus) or mostly (in ES) controlled by density, the $k_{\text{snow}}$ vertical profiles are also inverted, as already evidenced in Domine et al. (2016a). For example in early May 2015 (Fig. 4), we measured $k_{\text{snow}} = 0.028 \text{ W m}^{-1} \text{ K}^{-1}$ 2 cm above the ground, while Crocus (run *wind*) indicates $0.24 \text{ W m}^{-1} \text{ K}^{-1}$ and ES $0.35 \text{ W m}^{-1} \text{ K}^{-1}$. At 24 cm, values obtained from measurements for $k_{\text{snow}}$, Crocus and ES are respectively 0.28, 0.07 and $0.14 \text{ W m}^{-1} \text{ K}^{-1}$. Figure 5 summarizes the $k_{\text{snow}}$ measurements performed during field campaigns, as a function of the corresponding measured densities and snow types. Simulated values obtained on 14 May 2014 and 6 May 2015 from runs *wind* and *ES* are also shown. It confirms that high den-

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Measured and simulated snow thermal conductivity values as a function of density. Crocus computes $k_{\text{snow}}$ from the density only following Yen's parameterization, while ES includes the thermal effects of latent heat fluxes within the snowpack. Regression curves from Yen (1981) and Sturm et al. (1997) are also shown as well as observed snow types.</figcaption></figure>

sity values (> 400 kg m⁻³), mostly observed in wind slab layers, are neither reproduced by Crocus nor ES. For lower densities, it is well visible that density-$k_{\text{snow}}$ correlations (Eqs. 1 and 2) are not appropriate for our dataset. The simulated $k_{\text{snow}}$ values are almost always higher than measurements, and accounting for the Sun et al. (1999) additional term (run ES) amplifies the error. The regression curve from Sturm et al. (1997) fits our measurements better, because this parameterization is based on Arctic and subarctic snows instead of focusing on alpine conditions. However, for a given density, the $k_{\text{snow}}$ value can vary by a factor of 4 to 5 (Domine et al., 2016a; Sturm et al., 1997) so that density-$k_{\text{snow}}$ correlations cannot be used to accurately determine the thermal conductivity of Arctic snow. In fact, as demonstrated by Löwe et al. (2013), snow thermal conductivity depends on both density and microstructure. By adding second-order bounds based on a microstructural anisotropy parameter, they could improve the estimation of $k_{\text{snow}}$. In addition, our Fig. 5 strikingly corresponds to the theoretical Fig. 4 of Calonne et al. (2014). Our depth hoar values correspond to their lower bound, and our rounded grain snows (including wind slabs) correspond to their upper bound. Our data therefore confirms theoretical considerations which clearly demonstrate that parameterizing thermal conductivity as a function of density only simply should not be done, as it can produce very large errors.

### 3.1.3 Snow temperature

Snow density and $k_{\text{snow}}$ are the variables controlling the heat transfer through the snowpack, but they are significantly er-
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Evolution of simulated and observed temperatures within the snowpack during winter 2014–2015. Observations at 2, 12 and 22 cm come from the NPs, and the surface temperature is measured by the IR120 sensor.</figcaption></figure>

roneous in simulations. To evaluate the consequences on the temperature profile, Fig. 6 presents the evolution of temperature measured by the NPs at 2, 12 and 22 cm, the surface temperature measured by the infrared IR120 sensor and the corresponding simulated temperatures. The NPs perform measurements at 05:00 LT every other morning and record temperature only then. For this reason, simulated temperatures are shown every 2 days at about the same time. Because of the models' output resolution (6 h), simulated temperatures are actually shown at 07:00 LT, or 12:00 UTC.

The high variability in snow heights makes it difficult to estimate the actual snow accumulation around the NPs. But it appears that temperatures in the middle of the snowpack are reproduced best. This is well visible at 22 cm, with the temperature being very accurately simulated after January 2015, while the NP is definitively buried in snow. At the snow surface, the temperature is colder in simulations. This can be understood by considering that under a steady state approximation, the heat flux is the same through all snow layers and according to Fourier's law this flux is the product of $k_{\text{snow}}$ by the temperature gradient. Since above 22 cm the simulated $k_{\text{snow}}$ is lower than observed, the temperature gradient between 22 cm and the surface is greater and therefore the simulated surface temperature is colder. Following the same reasoning, the high simulated $k_{\text{snow}}$ at the bottom of the snowpack should lead to colder snow-soil interface temperature. However, temperatures simulated by Crocus are about 5 to 10 °C warmer than measured ones at the bottom of the snowpack. This is probably due to the fact that the overall thermal

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Overview of snow height (a), mean thermal conductivity (b) and thermal insulation (c) measured during May 2014 and 2015 field campaigns along with simulated values.</figcaption></figure>

conductivity of the snowpack is greater in Crocus simulations than observed (Fig. 5). The greater heat flux from the ground can then lead to a higher temperature at the base of the snow. For the ES simulation, the overall snow conductivity is so much higher than observed that ground cooling is much greater, as detailed below, and this leads to a colder temperature at the snow-ground interface.

An extended warm spell started on 16 March 2015, with the air temperature reaching –2.2 °C on 17 March associated with a significant snowfall (Fig. 2). This warmed up the entire snowpack, and from that date until the onset of simulated snowmelt, snow temperatures are well reproduced by the run *wind*. Snowmelt starts on 18 May, and as already discussed the snowpack is warming faster in simulations.

### 3.1.4 Sensitivity analysis on the snow thermal insulation

Properties of the snow pits studied in May 2014 and 2015 are summarized in Fig. 7, and compared to simulation outputs. In 2015, we used simulations of 6 May before the early melting. Since we are investigating the link between the soil temperature and snow thermal properties, Fig. 7 shows the snow height, the mean snow thermal conductivity $\bar{k}_{\text{snow}}$ and the resulting thermal insulation $R_T$ as given by Eqs. (3) and (4). Given that the NP method induces a systematic error, which underestimates $k_{\text{snow}}$ by about 20% on average, we also increased the measured $k_{\text{snow}}$ values by 20% (modified $k_{\text{snow}}$).
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Evolution of the soil thermal conductivity measured at 10 cm depth, and simulated (run *wind*) at depths of 10 and 20 cm.</figcaption></figure>

Litter and SOC additions have little effect on snow properties, the most noticeable being a reduction of less than 1 cm in snow height. Accounting for higher densities during wind compaction affects both the height and $k_{\text{snow}}$. Overall, the snow height is well reproduced in simulations and results are within the measured standard deviation obtained from snow pits. $\overline{k}_{\text{snow}}$ is overestimated for both years of simulations, resulting in a simulated thermal insulation lower than measured. Compared to measurements, the run *wind* gives a mean thermal conductivity 46% higher in 2014 and 73% higher in 2015, while the respective $R_T$ are 24 and 38% lower. The error in $R_T$ is thus essentially induced by the $k_{\text{snow}}$ simulation, and especially by the overestimation of the basal layer $k_{\text{snow}}$ value. Because ES simulates very high $k_{\text{snow}}$ values at the bottom of the snowpack (Fig. 4), the error in $R_T$ is amplified.

The larger error in $\overline{k}_{\text{snow}}$ simulated in the second winter can be explained by considering snow stratigraphies (Fig. 3). In May 2014, we observed a regular depth hoar basal layer whose density was successfully reproduced by Crocus. Its simulated $k_{\text{snow}}$ is still 4 times greater than measured (run *wind*). The following winter, the indurated depth hoar basal layer formed from a melted–refrozen snow. Crocus and ES do reproduce the partial melting of the snow at the beginning of the season, but cannot simulate the following transformation into depth hoar under high temperature gradients (Domine et al., 2016a). As the water vapor flux through the snowpack is not represented, the models still consider a refrozen layer at the bottom of the snowpack with a very high $k_{\text{snow}}$ (8 to 12 times greater than measured, Fig. 4). Hence, the mean simulated $k_{\text{snow}}$ is more overestimated in May 2015 than in 2014.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Vertical profiles of $k_{\text{soil}}$, soil temperature and volumetric water content (VWC, labeled Moisture in the figure) measured on 29 June 2014 and simulated in the first 20 cm below the surface. Horizontal bars indicate the standard deviation of measurements. The runs *litter*, *SOC*, *wind* and *ES* simulated the same soil temperature profiles.</figcaption></figure>

## 3.2 Simulations of soil properties and comparison with field data

### 3.2.1 Soil thermal conductivity

The monitoring of the soil thermal conductivity ($k_{\text{soil}}$) at –10 cm shows a bimodal distribution between the frozen and thawed state (Fig. 8). In summer, $k_{\text{soil}}$ is around 0.73 W m$^{-1}$ K$^{-1}$, while it suddenly increases with freezing to 1.95 W m$^{-1}$ K$^{-1}$ on average (Domine et al., 2016a). As detailed above, the noise in the frozen soil data is due to the low power used to heat the NPs. ISBA simulated $k_{\text{soil}}$ values from the run *wind* are shown at depths of 10 and 20 cm. A bimodal distribution is also visible with sharp transitions between the thawed and frozen state. The mean $k_{\text{soil}}$ value at –10 cm is 0.36 W m$^{-1}$ K$^{-1}$ in thawed conditions, and 2.04 W m$^{-1}$ K$^{-1}$ in the frozen state. The winter value is close to the measured mean, but the summer simulated $k_{\text{soil}}$ is lower than observations. Because the transition between litter and mineral soil occurs around 10 cm below the surface, measurements probably give a mixture of litter and mineral $k_{\text{soil}}$ values, while the simulated value at –10 cm is that of the litter only. Simulated $k_{\text{soil}}$ at 20 cm depth is 1.27 W m$^{-1}$ K$^{-1}$ in summer, and 1.89 W m$^{-1}$ K$^{-1}$ in winter. Averaging values between –10 and –20 cm results in a summer mean of 0.82 W m$^{-1}$ K$^{-1}$ and 1.97 W m$^{-1}$ K$^{-1}$ in winter, in very good agreement with our measurements. $k_{\text{soil}}$ is also greatly dependent on the water content, so we need to look at the vertical profiles of soil properties to assess the stratification.

Figure 9 shows vertical profiles of soil properties ($k_{\text{soil}}$, temperature and VWC) averaged from two measured profiles
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Observed and simulated daily mean soil temperature at 5, 10 and 15 cm deep.</figcaption></figure>

and simulated on 29 June 2014 to a depth of 20 cm. It illustrates the improvements caused by the litter addition in simulation: instead of staying constant at around 1.36 W m⁻¹ K⁻¹ in the absence of litter, $k_{\text{soil}}$ values are highly reduced in the first 10 cm when a litter is added, in fairly good agreement with our observations. However, we measured an increase in $k_{\text{soil}}$ with depth in the first 10 cm, from 0.19 W m⁻¹ K⁻¹ at -3 cm to 0.71 W m⁻¹ K⁻¹ at -10 cm, while the simulated $k_{\text{soil}}$ is constant (0.35 W m⁻¹ K⁻¹) through the litter layer. The same pattern is visible on the VWC profiles, attesting the water content dependence of $k_{\text{soil}}$. The lowest VWC value was found at -3 cm with 20% moisture, immediately increasing to more than 45% at -6 and -10 cm, while the simulated water content stays constant at 40% (runs *wind* and *ES*). The VWC difference at -10 cm is not sufficient to explain the $k_{\text{soil}}$ underestimation in summer, meaning that measurements of $k_{\text{soil}}$ are probably affected by both the litter and the underlying mineral soil. The water content starts to decrease below -10 cm, which supports the location of the lower limit of the litter. The presence of a litter in simulations finally improves the soil temperature profile between -10 and -20 cm, but the temperature is still too warm near the surface.

### 3.2.2 Soil temperature

Results of soil temperature simulations are shown in Fig. 10, along with measurement data from the Decagon sensors at depths of 5, 10 and 15 cm.

The run *base* simulates too high temperatures in both summer and winter, at all 3 depths shown. By reducing the heat exchanges between the atmosphere and the soil, litter and SOC additions greatly improve simulations in summer. Con-

sequences are well visible at -15 cm, with a good simulation of the observed temperature (+0.8 °C bias) during summer 2014 (June–July–August). For the same period, the temperature difference is +1.8 °C at -5 cm, and +1 °C at -10 cm. The same pattern is visible in Fig. 9, with a warm bias increasing toward the surface. The run *wind* improves the simulation of winter soil temperatures, resulting in the best simulation for the 2 years of measurements. Temperatures in winter months (December–January–February) are very well reproduced during 2014–2015, while during 2013–2014 there is a cold bias of 2.2 °C in December, which becomes less than 1 °C until snowmelt. ES produces soil temperatures up to 8 °C colder in winter, because it highly underestimates the snow thermal insulation (Fig. 7).

During the freezing and thawing periods, the models are not able to reproduce the observed soil temperature. The freezing of the active layer occurs in September and starts at the same time in simulations and observations, but it is too fast in simulations. The zero-curtain period lasts between 4 and 6 weeks in observations, but only a few days in simulations. This may be attributed to the thermal conductivity of the basal snow layer, which is too high in simulations, allowing the soil to cool rapidly. The thawing at the end of the winter is also faster in the simulations, again most likely because of the high thermal conductivity of the snow. The formation of a highly conductive refrozen snow layer is simulated at the top of the snowpack ($k_{\text{snow}}$ values up to 0.5 W m⁻¹ K⁻¹ are simulated in mid-May 2015 at 25 cm), while we measured values between 0.25 and 0.3 W m⁻¹ K⁻¹. Ultimately, the treatment of snow albedo in models is such that melt-out is hastened at high latitudes, inevitably accelerating soil thawing in spring.

Table 2 summarizes the performance of each run to reproduce the measured soil temperature at 15 cm depth. Depths of 5 and 10 cm give similar results. $R^2$ values are greater than 0.9, attesting the capacity of models to reproduce the observed variability. Statistics show a clear improvement as complexity is increased in the first 4 simulations. The presence of a litter noticeably improves the annual soil temperature simulation, while accounting for snow densities higher than 350 kg m⁻³ (run *wind*) drastically reduced the error in winter 2014–2015. The best results are obtained with the run *SOC* during 2013–2014 (+0.3 °C for the whole year), and with the run *wind* during 2014–2015 (+0.5 °C). Errors are very low for both years, of the order of the sensor accuracy.
<table><caption>**Table 2.** Statistical indicators of the differences between observations and model results of the soil temperature (°C) at -15 cm on a 6 h time step. Bold values correspond to the best scores.</caption><thead><tr><th rowspan="2">Year</th><th rowspan="2">period</th><th rowspan="2">type</th><th colspan="5">Runs</th></tr><tr><th>base</th><th>litter</th><th>SOC</th><th>wind</th><th>ES</th></tr></thead><tbody><tr><td rowspan="6">2013–2014</td><td rowspan="3">annual</td><td>bias</td><td>1.48</td><td>0.37</td><td><strong>0.32</strong></td><td>-0.92</td><td>-4.42</td></tr><tr><td>RMSE</td><td>3.41</td><td>2.62</td><td>2.61</td><td><strong>2.54</strong></td><td>5.91</td></tr><tr><td>r<sup>2</sup></td><td>0.92</td><td>0.94</td><td>0.94</td><td><strong>0.95</strong></td><td><strong>0.95</strong></td></tr><tr><td rowspan="3">winter (DJF)</td><td>bias</td><td>1.85</td><td>1.54</td><td><strong>1.33</strong></td><td>-1.42</td><td>-8.01</td></tr><tr><td>RMSE</td><td>1.94</td><td>1.66</td><td><strong>1.45</strong></td><td>1.62</td><td>8.06</td></tr><tr><td>r<sup>2</sup></td><td><strong>0.94</strong></td><td><strong>0.94</strong></td><td><strong>0.94</strong></td><td>0.88</td><td>0.93</td></tr><tr><td rowspan="6">2014–2015</td><td rowspan="3">annual</td><td>bias</td><td>3.31</td><td>2.08</td><td>2.04</td><td><strong>0.50</strong></td><td>-2.96</td></tr><tr><td>RMSE</td><td>4.34</td><td>3.49</td><td>3.52</td><td><strong>2.63</strong></td><td>4.86</td></tr><tr><td>r<sup>2</sup></td><td>0.94</td><td>0.94</td><td>0.94</td><td><strong>0.95</strong></td><td>0.92</td></tr><tr><td rowspan="3">winter (DJF)</td><td>bias</td><td>3.64</td><td>3.06</td><td>2.89</td><td><strong>0.38</strong></td><td>-5.07</td></tr><tr><td>RMSE</td><td>3.68</td><td>3.12</td><td>2.95</td><td><strong>0.61</strong></td><td>5.39</td></tr><tr><td>r<sup>2</sup></td><td>0.97</td><td>0.97</td><td>0.97</td><td><strong>0.98</strong></td><td>0.96</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Observed and simulated daily mean soil volumetric water content (VWC) at depths of 5, 10 and 15 cm.</figcaption></figure>

### 3.2.3 Soil water content

Figure 11 shows simulations of the soil VWC evolution at 5, 10 and 15 cm depths. The main difference with observation is the timing of freezing and thawing. This is clearly influenced by temperature, so that the errors in simulated temperature impact this timing. However, the duration of latent heat exchanges are also determined by the water content.

The VWC simulation in summer are improved by the addition of litter and SOC, but are still too low at -5 and -10 cm. Disabling the surface runoff helped to conserve high water contents in summer, but the moisture peak observed on 31

July 2014 at -5 cm is not simulated. In the ERAi precipitation data, a rain event did occur that day but it was less than 3 mm h⁻¹, not sufficient to increase the VWC by 20% as observed. Because of the high uncertainty on ERAi precipitation data, the actual precipitation rate could be underestimated. Figure 9 shows that the simulated VWC is constant in the first 10 cm of soil, while the observed profile is more variable. Thus, it seems that the simulation of the water dynamics in the first centimeter of the soil is erroneous, and is better reproduced below -10 cm. The low simulated VWC values therefore partly explain the excessively fast simulated soil freezing. However, as discussed by Langer et al. (2013), the inadequate thermal properties of the snow cover is most likely the main source of error in the ground thermal regime.

In winter, the different runs have a low impact on the water content, which is very well reproduced with around 6% of water remaining liquid while the temperature reaches -30 °C. This effect is explained by surface tensions applied on water confined in small pores, lowering its freezing temperature (Penner, 1970). The 6% threshold observed here is consistent with values found in the literature for a mixture of silt and fine sand.

Table 3 presents errors for each run in the VWC simulation, at the depth of 15 cm. Results at depths of 5 and 10 cm are similar. R² are lower than those from temperature simulations, ranging from 0.39 for the whole year (run base) to more than 0.9 in winter. SOC addition greatly improves simulations; it results in the lowest errors for the two winters and the lowest bias for annual simulations. Errors in winter are lower than the sensor resolution (3%), while annual simulations give RMSE greater than 7%, reflecting the limit of the model to reproduce the summer moisture variability.
<table><caption>**Table 3.** Statistical indicators of the differences between observations and model results of the soil volumetric water content (%) at –15 cm on a 6 h time step. Bold values indicate the best scores.</caption><thead><tr><th rowspan="2">Year</th><th rowspan="2">period</th><th rowspan="2">type</th><th colspan="5">Runs</th></tr><tr><th>base</th><th>litter</th><th>SOC</th><th>wind</th><th>ES</th></tr></thead><tbody><tr><td rowspan="6">2013–2014</td><td rowspan="3">annual</td><td>bias</td><td>-3.18</td><td>-4.18</td><td><strong>-2.21</strong></td><td>-2.38</td><td>-3.20</td></tr><tr><td>RMSE</td><td>11.85</td><td>8.90</td><td>8.37</td><td>8.19</td><td><strong>8.01</strong></td></tr><tr><td>r<sup>2</sup></td><td>0.48</td><td>0.73</td><td>0.72</td><td>0.74</td><td><strong>0.77</strong></td></tr><tr><td rowspan="3">winter (DJF)</td><td>bias</td><td>-2.97</td><td>-2.99</td><td><strong>-1.67</strong></td><td>-1.83</td><td>-2.13</td></tr><tr><td>RMSE</td><td>2.97</td><td>2.99</td><td><strong>1.67</strong></td><td>1.83</td><td>2.14</td></tr><tr><td>r<sup>2</sup></td><td><strong>0.93</strong></td><td><strong>0.93</strong></td><td><strong>0.93</strong></td><td>0.85</td><td>0.90</td></tr><tr><td rowspan="6">2014–2015</td><td rowspan="3">annual</td><td>bias</td><td>-2.31</td><td>-3.66</td><td><strong>-1.41</strong></td><td>-1.97</td><td>-2.69</td></tr><tr><td>RMSE</td><td>12.30</td><td>8.33</td><td>7.97</td><td><strong>7.29</strong></td><td>7.99</td></tr><tr><td>r<sup>2</sup></td><td>0.39</td><td>0.75</td><td>0.72</td><td><strong>0.78</strong></td><td>0.74</td></tr><tr><td rowspan="3">winter (DJF)</td><td>bias</td><td>-1.81</td><td>-1.85</td><td><strong>-0.52</strong></td><td>-0.67</td><td>-0.93</td></tr><tr><td>RMSE</td><td>1.81</td><td>1.85</td><td><strong>0.53</strong></td><td>0.67</td><td>0.93</td></tr><tr><td>r<sup>2</sup></td><td>0.89</td><td><strong>0.90</strong></td><td><strong>0.90</strong></td><td><strong>0.90</strong></td><td>0.89</td></tr></tbody></table>

# 4 Discussion

## 4.1 Snowpack simulations

In the Arctic, snowpacks are subjected to upward water vapor fluxes generated by strong temperature gradients between the atmosphere and the ground. These fluxes lead to mass transfers from the lower (warmer) to the upper (colder) snow layers. Consequently, this process has a considerable contribution to snow metamorphism by decreasing densities at the base of the snowpack and increasing densities of its upper parts (Domine et al., 2016a; Sturm and Benson, 1997). Domine et al. (2016a) estimated vertical water vapor fluxes in the snowpack, and came up with a mass loss of the basal layer of 2.6 kg m⁻² over two months. It corresponds to a density decrease from 300 to 200 kg m⁻³ for a 3 cm thick layer, in line with the model overestimation of the basal layer density. Although this estimation is approximate, it supports the suggestion that the vertical water vapor flux is the main cause of the model misrepresentation of the density profile. The resulting vertical $k_{\text{snow}}$ profiles are inverted in simulations, with high $k_{\text{snow}}$ at the bottom part of the snowpack and low $k_{\text{snow}}$ in the upper section, which affects the temperature gradient in the snowpack and the boundary fluxes. Further, downwelling shortwave absorption is lower in the Arctic because of the large zenith angle in late winter, and this is not accounted for in the original version of Crocus. Therefore, solar warming of the snowpack is exaggerated, resulting in incorrectly simulated melting episodes.

Despite their significant variability, averaged snow heights obtained from snow pits in May 2014 and 2015 are well reproduced. To obtain these, the snowfall amount in winter 2013–2014 had to be reduced by 30%. Without this artificial modification, the models were not able to reproduce this low

snow year. Given the importance of blowing snow events in the Arctic, accounting for snow redistribution by wind could also improve the simulated snow height (Gisnas et al., 2016; Libois et al., 2014). Even if relatively few strong wind events were recorded on Bylot Island in winter, Sturm et al. (2001b) showed that a single event could transport important amounts of snow. Hence, snow compaction and sublimation caused by wind are also critical to accurately simulate snow height and density in the Arctic, which was already improved after increasing the maximum snow density reached by drifting snow to 600 kg m⁻³. This change also partially compensates for vertical water transport, which increases upper layer densities. However, the models are not able to reproduce wind slab formation and the resulting increase in density as observed. But as already mentioned by Domine et al. (2016a): “The Crocus representation of the wind-packing process cannot be evaluated here, as the density increase also has contributions from water vapor deposition due to the upward flux and their respective contributions cannot be observed separately.”

## 4.2 Soil simulations

Granulometry is the main factor influencing soil thermal and hydraulic properties, so it is essential to analyze the soil composition of the studied site. Including SOC and especially litter greatly improved simulations of soil temperature and water content. Summer soil temperature is better reproduced at 15 cm depth than close to the surface. The addition of a litter improves simulation results, but the model is still restrained by a limited representation of surface organic covers. A detailed analysis of measured and simulated vertical profiles of soil properties highlighted that the difference between the surface moss and the dead litter leads to a stratification in the
first 2–3 cm of soil, with a lower thermal conductivity (i.e., a higher insulating capacity) for the moss. A moss surface also enhances the water infiltration (Beringer et al., 2001), resulting in a dryer soil surface and a greater storage of moisture in lower layers. Because $k_{\text{soil}}$ and the freezing process are highly dependent on the water content, improving the hydraulic properties of the soil scheme by considering a moss cover could help better reproduce the zero-curtain period and the soil thermal regime for layers close to the surface.

In light of the results of the snow–soil coupled evolution, it appears that errors in simulating the soil thermal regime manifest themselves mostly during freezing and thawing. In a similar study, Langer et al. (2013) found that the permafrost thermal state was mostly governed by the snow thermal properties. Here, the excessively rapid simulated soil freezing and thawing may be caused by the high thermal conductivity of the basal snow layer. For thawing, this is enhanced because of the early simulated snowmelt, linked to inadequate albedo simulations, which have a critical impact in May. It is interesting that the most sophisticated model experiment (run *wind*) is able to reproduce quite accurately the ground thermal regime at all depths in winter (Fig. 10), even though the thermal properties of the snow are not accurately simulated. Numerical models can be viewed as descriptions of a set of complex processes where error compensation is optimized. Therefore, we suggest that the insufficient insulating properties of the simulated snowpack are compensated by the inverted thermal conductivity profile, because the simulated insulating top snow layer damps air temperature fluctuations and thus limits heat transfer during cold spells. Under periodic diurnal temperature fluctuations, cold waves cannot penetrate into the simulated snowpack as effectively as when a conducting layer is present, so that the overall heat loss would be reduced by the inverted stratification. Of course, this is just a hypothesis that needs to be tested in future work, but some process exists, that compensates for the insufficiently insulating simulated snowpack, to explain the excellent soil temperature simulation most of the time.

## 4.3 Representativity of observations

The very large spatial variability of snow properties makes the comparison between simulation outputs and observations difficult. In particular, the snow height appeared to vary a lot with the microtopography on the 50 cm scale, thus only measurements performed in the immediate vicinity of the station can be used to assess the link between snow and soil properties measured at one specific site.

The lack of precipitation measurements is a major concern, because most of the snow and soil properties depend on the amount of precipitation. In particular, precipitation controls the snow height and the soil water content. It is well known that measuring precipitation in the Arctic is a challenge, especially in winter when snow falls in windy conditions and precipitation gauges are only catching between

20 and 70% of the actual amount of snow (Goodison et al., 1998; Liston and Sturm, 2002). But using a shielded precipitation gauge and correcting data using wind speed data should help reduce the large uncertainty that we have using ERA-Interim precipitation data alone, as already recommended by Bokhorst et al. (2016).

## 5 Conclusion

Applying the coupled snowpack-soil models ISBA-Crocus and ISBA-ES to a high Arctic site reveals major deficiencies related to typical Arctic conditions. The main weakness lies in the simulation of snow physical properties, because the absence of a modeled upward water vapor flux prevents reproducing the observed density profiles. The resulting vertical profiles of $k_{\text{snow}}$ are inverted in simulations, producing erroneous heat transfers through the snowpack. This work also illustrates that determining snow thermal conductivity from density only is inadequate especially (but not only) in the Arctic, confirming the theoretical work of Löwe et al. (2013) and Calonne et al. (2014). Considering also a microstructural variable, or perhaps snow type, appears mandatory (Lehning et al., 2002).

The soil temperature is the least well simulated during freezing and thawing periods. The main reason is the excessively conductive basal snow layer, which allows the soil to cool and warm too rapidly. Still, ISBA-Crocus manages to reproduce the temperature of the soil satisfactorily in summer and winter. Our results suggest that errors in the $k_{\text{snow}}$ stratification can compensate errors in simulated $k_{\text{snow}}$ values in winter. Hence, despite its apparent good results, Crocus is not better adapted than simpler models like ES to simulate Arctic snow thermal properties. The snow height has also a major influence on the winter soil temperature, but it is highly variable because of the wind-induced snow redistribution. Finally, a better representation of surface organic layers should improve simulations of the top soil properties, in particular the water content which controls the soil thermal properties and the water phase changes. The water content also governs the amount of water vapor transferred from the soil to the snow.

There are therefore strong uncertainties in climate projections related to the permafrost–carbon feedback, because even the most sophisticated snow models cannot accurately simulate heat transfers through Arctic snowpacks. In particular, the soil freezing and thawing processes need to be carefully simulated, because they determine most of the permafrost properties (Hinzman et al., 1991). The dual challenge to global climate models is thus to improve the representation of snow Arctic processes, in particular the water vapor fluxes and the resulting density and thermal conductivity profiles, and to include sub-grid snow height heterogeneities (Liston, 2004) in order to accurately simulate the thermal regime of permafrost.
*Code and data availability.* Forcing data used for this research and evaluation data of Figs. 2, 6, 8, 10 and 11 are available for further use in the scientific community at Nordicana D, https://doi.org/10.5885/45460CE-9B80A99D55F94D95 (Barrere and Domine, 2017). The strong differences between Arctic and alpine snowpack types encourage to incorporate such datasets in the upcoming model intercomparison exercises, in particular ESM-SnowMIP (http://www.geos.ed.ac.uk/~ressery/ESM-SnowMIP.html) which currently does not include sites featuring tundra snow types. The model is open source and it is available via the SURFEX platform which can be downloaded at http://www.umr-cnrm.fr/surfex/. The simulations presented in this paper were realized with SURFEX version 8.0, revision 4006.

*Author contributions.* FD and SM designed the research. MB and FD performed the field measurements and analyzed the data. MB realized the model simulations with advice from BD, SM, VV and ML. MB prepared the manuscript with contributions from all co-authors.

*Competing interests.* The authors declare that they have no conflict of interest.

*Acknowledgements.* This work was supported by the French Polar Institute (IPEV) through grant 1042 to Florent Domine, by NSERC through the Discovery Grants Program and by the BNP Paribas Foundation. Mathieu Barrere received scholarships from EnvironNord, and grants from Consulat Général de France à Québec and Région Rhône-Alpes. The Polar Continental Shelf Program (PCSP) efficiently provided logistical support for the research at Bylot Island. We are grateful to Denis Sarrazin for his precious technical support on the CEN weather stations, and to Gilles Gauthier and Marie-Christine Cadieux for their decades-long efforts to build and maintain the research base of the Centre d'Etudes Nordiques at Bylot Island. Bylot Island is located within Sirmilik National Park, and we thank Parks Canada and the Pond Inlet community (Mittimatilik) for permission to work there. CNRM/CEN and IGE are part of LabEx OSUG@2020. We thank the two anonymous reviewers for useful feedback and suggestions.

Edited by: David Lawrence

Reviewed by: two anonymous referees

## References

* ADAPT 2014: Carbon, nitrogen and water content of the active layer from sites across the Canadian Arctic, v. 1.0, Nordicana D20, https://doi.org/10.5885/45327AD-5245D08606AB4F52, 2014.

* Allard, M.: Geomorphological changes and permafrost dynamics: Key factors in changing Arctic ecosystems. An example from Bylot Island, Nunavut, Canada, Geosci. Canada, 23, 205–224, 1996.

* Barrere, M. and Domine, F.: Snow, soil and meteorological data at Bylot Island for simulating the permafrost thermal regime and evaluating output of the SURFEXv8 land surface scheme, v. 1.0 (1979–2015), Nordicana D29, https://doi.org/10.5885/45460CE-9B80A99D55F94D95, 2017.

* Bartelt, P. and Lehning, M.: A physical SNOWPACK model for the Swiss avalanche warning: Part I: numerical model, Cold Reg. Sci. Technol., 35, 123–145, https://doi.org/10.1016/S0165-232X(02)00074-5, 2002.

* Beringer, J., Lynch, A. H., Chapin, F. S. I., Mack, M., and Bonan, G. B.: The representation of arctic soils in the land surface model: The importance of mosses, J. Climate, 14, 3324–3335, https://doi.org/10.1175/1520-0442(2001)014<3324:TROASI>2.0.CO;2, 2001.

* Bokhorst, S., Pedersen, S. H., Brucker, L., Anisimov, O., Bjerke, J. W., Brown, R. D., Ehrich, D., Essery, R. L. H., Heilig, A., Ingvander, S., Johansson, C., Johansson, M., Jónsdóttir, I. S., Inga, N., Luojus, K., Macelloni, G., Mariash, H., McLenan, D., Rosqvist, G. N., Sato, A., Savela, H., Schneebeli, M., Sokolov, A., Sokratov, S. A., Terzago, S., Vikhamar-Schuler, D., Williamson, S., Qiu, Y., and Callaghan, T. V.: Changing Arctic snow cover: A review of recent developments and assessment of future needs for observations, modelling, and impacts, Ambio, 45, 516–537, https://doi.org/10.1007/s13280-016-0770-0, 2016.

* Boone, A., Masson, V., Meyers, T., and Noilhan, J.: The Influence of the Inclusion of Soil Freezing on Simulations by a Soil-Vegetation-Atmosphere Transfer Scheme, J. Appl. Meteorol., 39, 1544–1569, https://doi.org/10.1175/1520-0450(2000)039<1544:TIOTIO>2.0.CO;2, 2000.

* Boone, A. and Etchevers, P.: An Intercomparison of Three Snow Schemes of Varying Complexity Coupled to the Same Land Surface Model: Local-Scale Evaluation at an Alpine Site, J. Hydrometeorol., 2, 374–394, https://doi.org/10.1175/1525-7541(2001)002<0374:AIOTSS>2.0.CO;2, 2001.

* Brun, E., Martin, E., Simon, V., Gendre, C., and Coleou, C.: An energy and mass model of snow cover suitable for operational avalanche forecasting, J. Glaciol., 35, 333–342, 1989.

* Brun, E., Vionnet, V., Boone, A., Decharme, B., Peings, Y., Valette, R., Karbou, F., and Morin, S.: Simulation of Northern Eurasian Local Snow Depth, Mass, and Density Using a Detailed Snow-pack Model and Meteorological Reanalyses, J. Hydrometeorol., 14, 203–219, https://doi.org/10.1175/JHM-D-12-012.1, 2013.

* Calonne, N., Flin, F., Morin, S., Lesaffre, B., du Roscoat, S. R., and Geindreau, C.: Numerical and experimental investigations of the effective thermal conductivity of snow, Geophys. Res. Lett., 38, L23501, https://doi.org/10.1029/2011GL049234, 2011.

* Calonne, N., Flin, F., Geindreau, C., Lesaffre, B., and Rolland du Roscoat, S.: Study of a temperature gradient metamorphism of snow from 3-D images: time evolution of microstructures, physical properties and their associated anisotropy, The Cryosphere, 8, 2255–2274, https://doi.org/10.5194/tc-8-2255-2014, 2014.

* CEN 2016: Climate station data from Bylot Island in Nunavut, Canada, v. 1.7 (1992–2016), Nordicana D2, https://doi.org/10.5885/45039SL-EE76C1BDAADC4890, 2016.

* Charrois, L., Cosme, E., Dumont, M., Lafaysse, M., Morin, S., Libois, Q., and Picard, G.: On the assimilation of optical reflectances and snow depth observations into a detailed snowpack model, The Cryosphere, 10, 1021–1038, https://doi.org/10.5194/tc-10-1021-2016, 2016.

* Decharme, B., Boone, A., Delire, C., and Noilhan, J.: Local evaluation of the Interaction between Soil Biosphere At-
mosphere soil multilayer diffusion scheme using four pe-
dotransfer functions, J. Geophys. Res.-Atmos., 116, 1-29,
https://doi.org/10.1029/2011JD016002, 2011.

Decharme, B., Martin, E., and Faroux, S.: Reconciling soil thermal and hydrological lower boundary conditions in land surface models, J. Geophys. Res.-Atmos., 118, 7819–7834, https://doi.org/10.1002/jgrd.50631, 2013.

Decharme, B., Brun, E., Boone, A., Delire, C., Le Moigne, P., and Morin, S.: Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model, The Cryosphere, 10, 853–877, https://doi.org/10.5194/tc-10-853-2016, 2016.

Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. a., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597, https://doi.org/10.1002/qj.828, 2011.

Domine, F., Barrere, M., Sarrazin, D., Morin, S., and Arnaud, L.: Automatic monitoring of the effective thermal conductivity of snow in a low-Arctic shrub tundra, The Cryosphere, 9, 1265–1276, https://doi.org/10.5194/tc-9-1265-2015, 2015.

Domine, F., Barrere, M., and Sarrazin, D.: Seasonal evolution of the effective thermal conductivity of the snow and the soil in high Arctic herb tundra at Bylot Island, Canada, The Cryosphere, 10, 2573–2588, https://doi.org/10.5194/tc-10-2573-2016, 2016a.

Domine, F., Barrere, M., and Morin, S.: The growth of shrubs on high Arctic tundra at Bylot Island: impact on snow physical properties and permafrost thermal regime, Biogeosciences, 13, 6471–6486, https://doi.org/10.5194 bg-13-6471-2016, 2016b.

Elberling, B., Michelsen, A., Schädel, C., Schuur, E. A. G., Christiansen, H. H., Berg, L., Tamstorf, M. P., and Sigsgaard, C.: Long-term CO2 production following permafrost thaw, Nat. Clim. Change, 3, 890–894, https://doi.org/10.1038/nclimate1955, 2013.

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937–1958, https://doi.org/10.5194/gmd-9-1937-2016, 2016.

FAO/IIASA/ISRIC/ISSCAS/JRC: Harmonized World Soil Database (version 1.2), FAO, Rome, Italy and IIASA, Laxenburg, Austria, 2012.

Faroux, S., Le Moigne, P., and Dasprez, S.: http://www.umr-cnrm.fr/surfex/spip.php?rubrique143, last update: 3 May 2016.

Fortier, D. and Allard, M.: Late Holocene syngenetic ice-
wedge polygons development, Bylot Island, Canadian
Arctic Archipelago, Can. J. Earth Sci., 41, 997–1012,
https://doi.org/10.1139/e04-031, 2004.

Gisnas, K., Westermann, S., Vikhamar Schuler, T., Melvold, K., and Etzelmüller, B.: Small-scale variation of snow in a regional permafrost model, The Cryosphere, 10, 1201–1215, https://doi.org/10.5194/tc-10-1201-2016, 2016.

Goodison, B. E., Louie, P., and Yang, D.: WMO Solid Precipitation Measurement Intercomparison, WMO/TD, 872, 318 pp., 1998.

Gordon, M., Simon, K., and Taylor, P. A.: On snow depth predictions with the Canadian land surface scheme including a parametrization of blowing snow sublimation, Atmos. Ocean, 44, 239–255, https://doi.org/10.3137/ao.440303, 2006.

Guyomarc'h, G. and Mérindol, L.: Validation of an application for forecasting blowing snow, Ann. Glaciol., 26, 138–143, 1998.

Hinzman, L. D., Bettez, N. D., Bolton, W. R., Chapin, F. S., Dyurgierov, M. B., Fastie, C. L., Griffith, B., Hollister, R. D., Hope, A., Huntington, H. P., Jensen, A. M., Jia, G. J., Jorgenson, T., Kane, D. L., Klein, D. R., Kofinas, G., Lynch, A. H., Lloyd, A. H., McGuire, A. D., Nelson, F. E., Oechel, W. C., Osterkamp, T. E., Racine, C. H., Romanovsky, V. E., Stone, R. S., Stow, D. a., Sturm, M., Tweedie, C. E., Vourlitis, G. L., Walker, M. D., Walker, D. a., Webber, P. J., Welker, J. M., Winker, K. S., and Yoshikawa, K.: Evidence and Implications of Recent Climate Change in Northern Alaska and Other Arctic Regions, Clim. Change, 72, 251–298, https://doi.org/10.1007/s10584-005-5352-2, 2005.

Hinzman, L. D., Kane, D. L., Gieck, R. E., and Everett, K. R.: Hydrologic and thermal properties of the active layer in the Alaskan Arctic, Cold Reg. Sci. Technol., 19, 95–110, https://doi.org/10.1016/0165-232X(91)0001-W, 1991.

Hugelius, G., Strauss, J., Zubrzycki, S., Harden, J. W., Schuur, E. A. G., Ping, C.-L., Schirrmeister, L., Grosse, G., Michaelson, G. J., Koven, C. D., O'Donnell, J. A., Elberling, B., Mishra, U., Camill, P., Yu, Z., Palmtag, J., and Kuhry, P.: Estimated stocks of circumpolar permafrost carbon with quantified uncertainty ranges and identified data gaps, Biogeosciences, 11, 6573–6593, https://doi.org/10.5194/bg-11-6573-2014, 2014.

Jonasson, S., Chapin, F. S., and Shaver, G. R.: Biogeochemistry in the Arctic, Glob. Biogeochem. Cy.-Clim. Syst., 139–150, https://doi.org/10.1016/B978-012631260-7/50012-1, 2001.

Langer, M., Westermann, S., Heikenfeld, M., Dorn, W., and Boike, J.: Satellite-based modeling of permafrost temperatures in a tundra lowland landscape, Remote Sens. Environ., 135, 12–24, https://doi.org/10.1016/j.rse.2013.03.011, 2013.

Lawrence, D. M. and Slater, A. G.: The contribution of snow condi-
tion trends to future ground climate, Clim. Dynam., 34, 969–981,
https://doi.org/10.1007/s00382-009-0537-4, 2010.

Lehning, M., Bartelt, P., Brown, B., Fierz, C., and Satyawali, P.: A physical SNOWPACK model for the Swiss avalanche warning: Part II. Snow microstructure, Cold Reg. Sci. Technol., 35, 147–167, https://doi.org/10.1016/S0165-232X(02)00073-3, 2002.

Libois, Q., Picard, G., Arnaud, L., Dumont, M., Lafaysse, M., Morin, S., and Lefebvre, E.: Summertime evolution of snow specific surface area close to the surface on the Antarctic Plateau, The Cryosphere, 9, 2383–2398, https://doi.org/10.5194/tc-9-2383-2015, 2015.

Libois, Q., Picard, G., Arnaud, L., Morin, S., and Brun, E.: Modeling the impact of snow drift on the decameter-scale variability of snow properties on the Antarctic Plateau, J. Geophys. Res.-Atmos., 119, 11662–11681, https://doi.org/10.1002/2014JD022361, 2014.

Liston, G. E.: Representing Subgrid Snow Cover Heterogeneities in Regional and Global Models, J. Climate, 17, 1381–1397, https://doi.org/10.1175/1520-0442(2004)017<1381:RSSCHI>2.0.CO;2, 2004.
Liston, G. E. and Sturm, M.: Winter Precipitation Patterns in Arctic Alaska Determined from a Blowing-Snow Model and Snow-Depth Observations, J. Hydrometeorol., 3, 646–659, https://doi.org/10.1175/1525-7541(2002)003<0646:WPPIAA>2.0.CO;2, 2002.

Löwe, H., Riche, F., and Schneebeli, M.: A general treatment of snow microstructure exemplified by an improved relation for thermal conductivity, Cryosph., 1473–1480, https://doi.org/10.5194/tc-7-1473-2013, 2013.

Masson, V., Le Moigne, P., Martin, E., Faroux, S., Alias, A., Alkama, R., Belamari, S., Barbu, A., Boone, A., Bouyssel, F., Brousseau, P., Brun, E., Calvet, J.-C., Carrer, D., Decharme, B., Delire, C., Donier, S., Essaouini, K., Gibelin, A.-L., Giordani, H., Habets, F., Jidane, M., Kerdraon, G., Kourzeneva, E., Lafaysse, M., Lafont, S., Lebeaupin Brossier, C., Lemonsu, A., Mahfouf, J.-F., Marguinaud, P., Mokhtari, M., Morin, S., Pigeon, G., Salgado, R., Seity, Y., Taillefer, F., Tanguy, G., Tulet, P., Vincendon, B., Vionnet, V., and Voldoire, A.: The SURFEXv7.2 land and ocean surface platform for coupled or offline simulation of earth surface variables and fluxes, Geosci. Model Dev., 6, 929–960, https://doi.org/10.5194/gmd-6-929-2013, 2013.

Morin, S., Domine, F., Arnaud, L., and Picard, G.: In-situ monitoring of the time evolution of the effective thermal conductivity of snow, Cold Reg. Sci. Technol., 64, 73–80, https://doi.org/10.1016/j.coldregions.2010.02.008, 2010.

Noilhan, J. and Mahfouf, J. F.: The ISBA land surface parameterisation scheme, Global Planet. Change, 13, 145–159, https://doi.org/10.1016/0921-8181(95)00043-7, 1996.

Noilhan, J. and Planton, S.: A Simple Parameterization of Land Surface Processes for Meteorological Models, Mon. Weather Rev., 117, 536–549, https://doi.org/10.1175/1520-0493(1989)117<0536:ASPOLS>2.0.CO;2, 1989.

Pearson, R. G., Phillips, S. J., Loranty, M. M., Beck, P. S. A., Damoulas, T., Knight, S. J., and Goetz, S. J.: Shifts in Arctic vegetation and associated feedbacks under climate change, Nat. Clim. Change, 3, 673–677, https://doi.org/10.1038/nclimate1858, 2013.

Peters-Lidard, C. D., Blackburn, E., Liang, X., and Wood, E. F.: The effect of soil thermal conductivity parameterization on surface energy fluxes and temperatures, J. Atmos. Sci., 55, 1209–1224, 1998.

Post, E. and Pedersen, C.: Opposing plant community responses to warming with and without herbivores, P. Natl. Acad. Sci., 105, 12353–12358, https://doi.org/10.1073/pnas.0802421105, 2008.

Post, E., Forchhammer, M. C., Bret-Harte, M. S., Callaghan, T. V., Christensen, T. R., Elberling, B., Fox, A. D., Gilg, O., Hik, D. S., Høye, T. T., Ims, R. A., Jeppesen, E., Klein, D. R., Madsen, J., McGuire, A. D., Rysgaard, S., Schindler, D. E., Stirling, I., Tamstorf, M. P., Tyler, N. J. C., van der Wal, R., Welker, J., Wookey, P. A., Schmidt, N. M., and Aastrup, P.: Ecological dynamics across the Arctic associated with recent climate change, Science, 325, 1355–1358, https://doi.org/10.1126/science.1173113, 2009.

Riche, F. and Schneebeli, M.: Thermal conductivity of snow measured by three independent methods and anisotropy considerations, The Cryosphere, 7, 217–227, https://doi.org/10.5194/tc-7-217-2013, 2013.

Schaefer, K., Lantuit, H., Romanovsky, V. E., Schuur, E. G. and Witt, R.: The impact of the permafrost carbon feedback

on global climate, Environ. Res. Lett., 9, 085003, https://doi.org/10.1088/1748-9326/9/8/085003, 2014.

Schuur, E. A. G., McGuire, A. D., Grosse, G., Harden, J. W., Hayes, D. J., Hugelius, G., Koven, C. D., and Kuhry, P.: Climate change and the permafrost carbon feedback, Nature, 520, 171–179, https://doi.org/10.1038/nature14338, 2015.

Serreze, M. C., Walsh, J. E., Chapin, F. S. I., Osterkamp, T., Dyurgerov, M., Romanovsky, V., Oechel, W. C., Morison, J., Zhang, T., and Barry, R. G.: Observational evidence of recent change in the northern high latitude environment, Clim. Change, 46, 159–207, https://doi.org/10.1023/A:1005504031923, 2000.

Smith, S. L., Romanovsky, V. E., Lewkowicz, A. G., Burn, C. R., Allard, M., Clow, G. D., Yoshikawa, K., and Throop, J.: Thermal state of permafrost in North America: a contribution to the international polar year, Permafr. Periglac., 21, 117–135, https://doi.org/10.1002/ppp.690, 2010.

Sturm, M. and Benson, C. S.: Vapor transport, grain growth and depth-hoar development in the subarctic snow, J. Glaciol., 43, 42–59, 1997.

Sturm, M., Holmgren, J., König, M., and Morris, K.: The thermal conductivity of seasonal snow, J. Glaciol., 43, 26–41, 1997.

Sturm, M., McFadden, J. P., Liston, G. E., Stuart Chapin, F., Racine, C. H., and Holmgren, J.: Snow-shrub interactions in Arctic Tundra: A hypothesis with climatic implications, J. Climate, 14, 336–344, https://doi.org/10.1175/1520-0442(2001)014<0336:SSIIAT>2.0.CO;2, 2001a.

Sturm, M., Liston, G. E., Benson, C. S., and Holmgren, J.: Characteristics and Growth of a Snowdrift in Arctic Alaska, U.S.A., Arctic, Antarct. Alp. Res., 33, 319–329, https://doi.org/10.2307/1552239, 2001b.

Sun, S., Jin, J., and Xue, Y.: A simple snow-atmosphere-soil transfer model, J. Geophys. Res., 104, 19587, https://doi.org/10.1029/1999JD900305, 1999.

Sweet, S. K., Griffin, K. L., Steltzer, H., Gough, L., and Boelman, N. T.: Greater deciduous shrub abundance extends tundra peak season and increases modeled net CO₂ uptake, Global Change Biol., 21, 2394–2409, https://doi.org/10.1111/gcb.12852, 2015.

Tuzet, F., Dumont, M., Lafaysse, M., Picard, G., Arnaud, L., Voisin, D., Lejeune, Y., Charrois, L., Nabat, P., and Samuel, M.: A multi-layer physically-based snowpack model simulating direct and indirect radiative impacts of light-absorbing impurities in snow, The Cryosphere Discuss., https://doi.org/10.5194/tc-2017-94, in review, 2017.

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J. M.: The detailed snow-pack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791, https://doi.org/10.5194/gmd-5-773-2012, 2012.

Vionnet, V., Guyomarc'h, G., Naaim Bouvet, F., Martin, E., Durand, Y., Bellot, H., Bel, C., and Puglièse, P.: Occurrence of blowing snow events at an alpine site over a 10-year period: Observations and modelling, Adv. Water Resour., 55, 53–63, https://doi.org/10.1016/j.advwatres.2012.05.004, 2013.

Vuichard, N. and Papale, D.: Filling the gaps in meteorological continuous data measured at FLUXNET sites with ERA-Interim reanalysis, Earth Syst. Sci. Data, 7, 157–171, https://doi.org/10.5194/essd-7-157-2015, 2015.

Warren, S. G.: Optical properties of snow, Rev. Geophys., 20, 67–89, https://doi.org/10.1029/RG020i001p00067, 1982.
Yen, Y. C.: Review of the thermal properties of snow, ice and sea ice, Tech. Rep., 81-10, Cold Reg. Res. Eng. Lab., Hanover, NH, 1981.

Zhang, T.: Influence of seasonal snow cover on the ground thermal regime: an overview, Rev. Geophys., 43, RG4002, https://doi.org/10.1029/2004RG000157, 2005.