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

Global warming projections still suffer from a limited representation of the permafrost-carbon feedback. Predicting the response of the permafrost temperature to climate changes requires accurate simulations of the Arctic snow and soil 15 properties. This study assesses the capacity of the coupled 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 reanalysis were used to drive the models and to evaluate simulation outputs. Snow height, density, temperature, thermal conductivity and thermal resistance are examined to determine the critical variables involved in the soil thermal regime. Simulated soil properties are compared with measurements of thermal conductivity, temperature and water content. The simulated snow density profiles are 20 erroneous, because Crocus and ES do not represent the upward water vapour 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 representing adequately surface organic layers, i.e. mosses and litter, and in particular their water retention capacity. Transition 25 periods (soil freezing and thawing) are the least well reproduced because the high basal snow thermal conductivity induces too rapid heat transfers 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 changes. 30 Geosci. Model Dev. Discuss., doi:10.5194/gmd-2017-50, 2017 Manuscript under review for journal Geosci. Model Dev. Discussion started: 21 March 2017 c © Author(s) 2017. CC-BY 3.0 License.

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

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

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 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. (2015Domine et al. ( , 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 snow values (0.025-0.7 W m −1 K −1 ), 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 3 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 soil .The method used is the same as for the k 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).

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.

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 longwave 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 −4 g kg −1 in air humidity and +0.7 m s −1 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-similation 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.

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 icerich 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 −1 K −1 .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 −1 K −1 when dry, and to 0.25 W m −1 K −1 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 −3 for the first soil horizon (0-30 cm), and 4.34 kg m −3 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 −3 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).

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 −1 (Sturm et al., 2001b;Vionnet et al., 2013).From the 10 m tower on Bylot Island, wind speeds greater than 5 m s −1 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 onedimensional (Brun et al., 2013).In the standard version of Crocus, wind-compacted snow layers can reach a maximum density of 350 kg m −3 .As we observed snow densities up to 450 kg m −3 on Bylot Island, and values up to 600 kg m −3 can be found in the literature (e.g., Sturm et al., 1997;Zhang, 2005), we increased this value to 600 kg m −3 (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 snow , in W m −1 K −1 , which is used to solve the thermal diffusion equation in the snowpack, is calculated from the density using the equation of Yen (1981): with k ice the thermal conductivity of ice (2.22 W m −1 K −1 ), ρ the density of snow and ρ w the density of liquid water (1000 kg m −3 ).
www.geosci-model-dev.net/10/3461/2017/Geosci.Model Dev., 10, 3461-3479, 2017 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, lightabsorbing spectroscopy) prevented us from using TARTES on Bylot Island simulations.

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 snow calculation method is different from Crocus.In ES, k 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): 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 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.

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.

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 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 insulance of the snowpack R T (in m 2 K W −1 ) depends on the thickness h and the thermal conductivity of each layer i: The simulated R T is computed from the respective densityk 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 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: 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.

bias
where sim and obs refer to simulated and observed values respectively.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 mea- www.geosci-model-dev.net/10/3461/2017/Geosci.Model Dev., 10, 3461-3479, 2017 sured by our instruments, and because the large spatial variability 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 −1 ), 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 simulation results of 6 May, just before the partial melting.

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 bottom depth hoar layer (between 150 and 200 kg m −3 ), and high values for the upper wind slabs (exceeding 400 kg m −3 ).Crocus and ES simulate inverted profiles, with generally decreasing densities from the bottom to the top of the snowpack.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 reproduced 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 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 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 snow = 0.028 W m −1 K −1 2 cm above the ground, while Crocus (run wind) indicates 0.24 W m −1 K −1 and ES 0.35 W m −1 K −1 .At 24 cm, values obtained from measurements for k snow , Crocus and ES are respectively 0.28, 0.07 and 0.14 W m −1 K −1 .Figure 5   sity values (> 400 kg m −3 ), mostly observed in wind slab layers, are neither reproduced by Crocus nor ES.For lower densities, it is well visible that density-k snow correlations (Eqs. 1 and 2) are not appropriate for our dataset.The simulated k 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 snow value can vary by a factor of 4 to 5 (Domine et al., 2016a;Sturm et al., 1997) so that density-k 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 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.

Snow temperature
Snow density and k snow are the variables controlling the heat transfer through the snowpack, but they are significantly erwww.geosci-model-dev.net/10/3461/2017/Geosci.Model Dev., 10, 3461-3479, 2017 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 snow by the temperature gradient.Since above 22 cm the simulated k 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 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 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.

Sensitivity analysis on the snow thermal insulance
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 k snow and the resulting thermal insulance R T as given by Eqs. ( 3) and ( 4).Given that the NP method induces a systematic error, which underestimates k snow by about 20 % on average, we also increased the measured k snow values by 20 % (modified k snow ).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 snow .Overall, the snow height is well reproduced in simulations and results are within the measured standard deviation obtained from snow pits.k snow is overestimated for both years of simulations, resulting in a simulated thermal insulance 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 snow simulation, and especially by the overestimation of the basal layer k snow value.Because ES simulates very high k snow values at the bottom of the snowpack (Fig. 4), the error in R T is amplified.
The larger error in k 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 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 snow (8 to 12 times greater than measured, Fig. 4).Hence, the mean simulated k snow is more overestimated in May 2015 than in 2014.3.2 Simulations of soil properties and comparison with field data

Soil thermal conductivity
The monitoring of the soil thermal conductivity (k soil ) at −10 cm shows a bimodal distribution between the frozen and thawed state (Fig. 8).In summer, k 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 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 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 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 soil values, while the simulated value at −10 cm is that of the litter only.Simulated k 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 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 soil , temperature and VWC) averaged from two measured profiles 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 −1 K −1 in the absence of litter, k 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 soil with depth in the first 10 cm, from 0.19 W m −1 K −1 at −3 cm to 0.71 W m −1 K −1 at −10 cm, while the simulated k soil is constant (0.35 W m −1 K −1 ) through the litter layer.The same pattern is visible on the VWC profiles, attesting the water content dependence of k 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 the k soil underestimation in summer, meaning that measurements of k 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.

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 insulance (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 snow values up to 0.5 W m −1 K −1 are simulated in mid-May 2015 at 25 cm), while we measured values between 0.25 and 0.3 W m −1 K −1 .Ultimately, the treatment of snow albedo in models is such that meltout 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 −3 (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.

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 −1 , 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 2 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.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).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 −3 .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."

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

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

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 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 snow stratification can compensate errors in simulated k 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.

Figure 1 .
Figure 1.Location of the study site in the southwest plain of Bylot Island in the Canadian Arctic archipelago.

Figure 2 .
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.

Figure 2
Figure2shows 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 Table2ofDomine 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 mea-

Figure 3 .
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.

Figure 4 .
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.
summarizes the k 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 5 .
Figure 5. Measured and simulated snow thermal conductivity values as a function of density.Crocus computes k snow from the density only following Yen's parameterization, while ES includes the thermal effects of latent heat fluxes within the snowpack.Regression curves fromYen (1981) andSturm et al. (1997)  are also shown as well as observed snow types.

Figure 6 .
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.

Figure 7 .
Figure 7. Overview of snow height (a), mean thermal conductivity (b) and thermal insulance (c) measured during May 2014 and 2015 field campaigns along with simulated values.

Figure 8 .
Figure 8. Evolution of the soil thermal conductivity measured at 10 cm depth, and simulated (run wind) at depths of 10 and 20 cm.

Figure 9 .
Figure 9. Vertical profiles of k 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.

Figure 10 .
Figure 10.Observed and simulated daily mean soil temperature at 5, 10 and 15 cm deep.

Figure 11 .
Figure 11.Observed and simulated daily mean soil volumetric water content (VWC) at depths of 5, 10 and 15 cm.

Table 1 .
Configurations of the numerical experiments, all realized with Crocus except for ES.Changes are incremented in addition to the previous configuration.

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.

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.
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 −2 over two months.It corresponds to a density decrease from 300 to 200 kg m −3 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 snow profiles are inverted in simulations, with high k snow at the bottom part of the snowpack and low k 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.