Impacts of snow and organic soils parameterization on North-Eurasian soil temperature profiles

Introduction Conclusions References


Introduction
Snowpack properties are known to be of primary importance for understanding the water and energy budgets of the land surface, especially in mountainous and boreal regions. From autumn to spring, solid precipitation is stored within the snowpack thereby modifying the terrestrial albedo and roughness length, and impacting the this last process controls the temperature of the permafrost. It is defined as a soil that remains below 0 • C for two or more consecutive years, and it has a significant influence on the summer permafrost active layer thickness, defined as the maximum annual thaw depth. In summary, snowpack properties drastically influence soil/atmosphere interactions during a large part of the year through their impacts on many land surface 10 processes.
Beside the importance of snowpack properties for understanding the water and energy budgets of the land surface in northern regions, the physical properties of soil organic carbon (or peat soil) also play a significant role. North-Eurasian soils are very rich in organic carbon because the low soil temperatures in this region inhibit 15 decomposition of dead plant material that accumulates over time, thereby forming peat deposits. Soil organic carbon exhibits very different hydraulic and thermal properties than mineral soil (Boelter, 1969;Letts et al., 2000). It is characterized by a very high porosity, a weak hydraulic suction, and a sharp vertical hydraulic conductivity profile from high values at the surface to very low values at the subsurface. This generally 20 induces a relatively wet soil with a shallow water table (Letts et al., 2000). Its low thermal conductivity and its relatively high heat capacity act as an insulator for soil temperature that prevents the soil from significant warming during the summer (Bonan and Shugart, 1989;Lawrence and Slater, 2008). Over permafrost regions, the hydraulic and thermal properties of soil organic carbon partly control the soil depth reached by 25 the 0 • C isotherm which, in turn, defines the thickness of the active layer during summer (Paquin and Sushama, 2015). Through its influence on soil temperature and wetness, it impacts the continental part of the carbon cycle and the land surface CO 2 and CH 4 emissions to the atmosphere (Walter et al., 2006;Zimov et al., 2006).

6735
Now, LSMs are used in many applications such as hydrological and meteorological forecasts, global hydrological and biogeochemical studies, and climate evolution prediction. Many LSMs use multi-layer soil schemes in which the vertical transport of moisture and heat into the soil is explicitly solved for using diffusion equations (e.g. Decharme et al., 2011). Because the total soil depth is discretized using multiple layers, 10 these schemes allow the representation of the vertical root distribution (Zeng et al., 1998;Feddes et al., 2001;Braud et al., 2005), as well as the surface/groundwater capillary exchanges (e.g. Vergnes et al., 2014). Finally, their coupling with a multi-layer snowpack scheme permits a representation of the interaction between cold physical processes, such as the effect of snow on soil temperature, hydrology, and freezing 15 (Slater et al., 2001;Luo et al., 2003;Gouttevin et al., 2012). Three major classes of snowpack schemes exist in LSMs: single-layer schemes, multi-layer schemes of intermediate complexity, and detailed snowpack models. The first class was used preferentially in the past within forecast and climate models. The snowpack was represented with only one layer that evolves seasonally, which us systems, global hydrological models and model chains in support of avalanche hazard warning (e.g. Lafaysse et al., 2013;Vernay et al., 2015). The ISBA multi-layer version was evaluated over many local or regional field datasets (Boone et al., 2000;Decharme et al., 2011Decharme et al., , 2013Canal et al., 2014;Parrens et al., 2014;Vergnes et al., 2014;Joetzjer et al., 2015), increasing our confidence in the model's capability to simulate realistic 20 land surface processes under a variety of climate conditions. However, over cold regions, winter top soil temperatures tend to be underestimated (Wang et al., 2015) while during summer they are generally too warm. The first biases are attributable to the ISBA multi-layer snowpack scheme of intermediate complexity developed by Boone et al. (2000) and based on Anderson (1979). Indeed, when the ISBA multi-Introduction carbon are required to simulate realistic soil thermal regime over cold regions (Nicolsky et al., 2007;Beringer et al., 2001;Lawrence and Slater, 2008;Lawrence et al., 2008;Dankers et al., 2011). The present study focuses on the impact of improving the representation of snowpack and soil properties in the ISBA LSM to reproduce snow characteristics 5 and soil temperature profiles over cold regions. We replaced the original Boone and Etchevers (2001) representation of snow layering, albedo and snow compaction by adapting some parameterizations used in the Crocus snowpack model (e.g. Vionnet et al., 2012). In addition, we added a parameterization of the organic carbon effect on hydraulic and thermal soil properties based on the pedotransfer function of Boelter 10 (1969) and inspired by works of Letts et al. (2000) and Lawrence and Slater (2008). The changes in the snowpack parameterizations are first evaluated at the Col de Porte field site located in the French Alps . This dataset includes many observations at a daily time step such as snow depth, snow water equivalent, surface albedo and soil temperature at 10 cm from 1993 to 2011. In addition the meteorological 15 observations required to drive the model are given at a 3-hourly time step over the same period. The new parameterizations were evaluated next over the North-Eurasian region using the same experimental design as Brun et al. (2013) using in situ evaluation datasets of snow depth and soil temperature profile measurements and meteorological driving data from a global reanalysis. To quantify the model's ability to simulate the 20 permafrost characteristics, two additional datasets were used that estimate the location of permafrost boundaries and the active layer thickness over the Yakutia region. A brief review of the ISBA multi-layer model is given in Sect. 2, all of the snowpack and soil parameterization improvements and updates are presented in Sect. 3. Sections 4 and 5 describe the model evaluation over the Col de Porte field site and the North-Eurasian Introduction

Soil processes
The ISBA multi-layer model solves the one-dimensional Fourier law and the mixed-form of the Richards equation explicitly to calculate the time evolution of the soil energy and water budgets (Boone et al., 2000;Decharme et al., 2011). In each layer i , the 5 closed-form equations between the soil liquid water content, w (m 3 m −3 ), and the soil hydrodynamic parameters, such as the soil matric potential, ψ (m), and the hydraulic conductivity, k (m s −1 ), are determined according to the Brooks and Corey (1966) model adapted by Campbell (1974) as follows: (1) 10 where, b represents the dimensionless shape parameter of the soil-water retention curve, w sat (m 3 m −3 ) the soil porosity, and ψ sat (m) and k sat (m s −1 ) the soil matric potential and hydraulic conductivity at saturation, respectively. In this study, the heat and soil moisture transfers within the soil are computed using 14 layers up to a 12 m depth. The depth of the 14 layers (0.01, 0.04, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 15 3.0, 5.0, 8.0, 12.0 m) have been chosen to minimize numerical errors in solving the finite-differenced diffusive equations, especially in the uppermost meter of the soil . Saturated hydraulic conductivity, matric potential at saturation, and porosity of the mineral soil are related to the soil texture (Noilhan and Lacarrère 1995). The total heat capacity of the mineral soil in each layer is computed as the sum 20 of the soil matrix, water and ice heat capacities weighted by the volumetric water and ice content (Peters-Lidard et al., 1998). The thermal conductivity of the mineral soil is computed via a more complex combination of water, ice and soil conductivities as proposed by Peters-Lidard et al. (1998 surface (e.g. Boone et al., 2000). The liquid water content that can freeze is limited by a maximum value, w lmax (m 3 m −3 ), computed as a function of temperature based on the Gibbs free-energy method (Fuchs et al., 1978): where w sat (m 3 m −3 ) is the soil porosity in each layer i , T g (K) the soil temperature, g 5 (m s −2 ) the terrestrial gravity constant, T f (273.16 K) is the triple-point temperature for water, and L f (3.337 × 10 5 J kg −1 ) the latent heat of fusion. The total water content in each soil layer is conserved during phase changes. When the soil freezes, the liquid water content will decrease owing to a corresponding increase in soil ice content. Finally, the maximum temperature, T max (K), used for phase changes can be 10 determined via the same Gibbs free-energy method: where the soil matric potential ψ is defined using Eq. (1). Thus, this scheme induces dependencies of water phase changes to soil textures and to the degree of soil humidity. The coarser the soil texture, the larger the quantity of water that will freeze 15 at a given temperature. As the soil becomes dry, the temperature that allows freezing drops. More details can be found in the Supplement of Masson et al. (2013)  snow compaction, water percolation and refreezing, and explicit heat conduction at the snow/soil interface. Many of theses processes, such as snow compaction or absorption of solar energy, are based on works of Anderson (1976) and Loth et al. (1993). The thermal conductivity of snow (Appendix A) is computed via the snow density (Yen, 1981). An additional term depends on the snow temperature to account for vapor 5 transfer through the snowpack (Sun et al., 1999). The time evolution of the snow mass is linked to snowmelt, water freezing, evaporation, and liquid flow. The liquid water content into the snowpack is simulated as a succession of bucket-type reservoirs. A maximum liquid water holding capacity (W lmax ) is computed in each layer. It varies from 3 to 10 % of the snow mass according to a decrease in snow density after 10 Anderson (1976). A liquid water flux is generated when the liquid water content exceeds this threshold. More details can be found in Boone and Etchevers (2001) and only internal physical processes of the snowpack discussed in this study are described below. 15 In the original ISBA explicit snow scheme, three-layers are used for snow layering because it is considered to be the minimum number required to resolve adequately the snow thermal profile within the snowpack (Lynch-Stieglitz, 1994;Loth and Graf, 1998;Boone and Etchevers, 2001). The algorithm that computes the snow grid thicknesses ∆z of each layer, i , is described as follows: where h sn (m) is the total snow depth. As long as the snow remains below 0.2 m, the fraction of the total depth that defines the thickness of each layer remains with a fine resolution at the top and the base of the snowpack. When the snow depth exceeds 0.2 m, the thickness of the first layer remains equal to 0.05 m, in order to adequately solve the diurnal cycle of the surface energy balance. In addition, for large snow depth 5 values, the second layer thickness cannot exceed 0.5 m because density and heat vertical gradients are generally the largest near the top of the snowpack. The vertical grid is updated at the beginning of each time step before the computation of the other snowpack internal processes. 10 The evolution of snow density, ρ sn (kg m −3 ) in each layer, i , is the sum of snow compaction due to change in snow viscosity, η (Pa s −1 ), and settling due to freshly fallen snow, ξ (s −1 ), following Anderson (1976) and Loth et al. (1993):

Snow compaction
where σ (Pa) is the snow vertical stress. The snow viscosity and settling of new snow 15 are solved using two empirical exponential functions of snow density and temperature, T sn (K): where v 0 = 3.7×10 7 Pa s −1 , v 1 = 0.081 K −1 , v 2 = 0.018 m 3 kg −1 , s 0 = 2.8×10 −6 s −1 , s 1 = 20 0.04 K −1 , s 2 = 0.046 m 3 kg −1 , and ρ d = 150 kg m −3 are empirical parameters calibrated by Anderson (1976). The minimum density of snow is constrained to be 50 kg m −3 . The snowfall density, ρ snew (kg m −3 ), is expressed as a function of wind speed, V a (m s −1 ), and air temperature, T a (K), following an experimental study of Pahaut (1976): where the coefficients a ρ = 109 kg m −3 , b ρ = 6 kg m −3 K −1 , and c ρ = 26 kg s 1/2 m −7/2 . 5 The absorption of incident shortwave solar radiation, R SW (W m −2 ), within the snowpack is solved over a single spectral band. It uses an exponential decrease of incoming radiation with snow depth (Anderson, 1976;Loth et al., 1993). So, the net shortwave radiation Q sn (W m −2 ) absorbed by the snow level, i , is given by:

Transmission of solar radiation and snow albedo
where α sn is the dimensionless snow albedo, and β sn (m −1 ) the extinction coefficient of snow which is given by: As shown by Bohren and Barkstrom (1974), this extinction of snow is directly related to its density, the optical diameter d opt (m), and a constant C ν = 3.8 × 10 −3 m 5/2 kg −1 . The 15 optical diameter is empirically linked to the snow density following a simple polynomial regression established by Anderson (1976): where d max (m) is the maximum value equal to 2.796×10 −3 m, and the coefficients g 1 = 1.6 × 10 −4 m, and g 2 = 1.1 × 10 −13 m 13 kg −4 were calibrated by Anderson (1976). The time evolution of snow albedo is modelled in a simple way using time constants after Douville et al. (1995). A linear decrease rate is used for dry snow and an exponential decrease is used for wet snow. The snow albedo is constrained to be between its 5 minimum value, α min = 0.5, and its maximum, α max = 0.85.

Snow layering
Detailed snowpack models use more than a dozen layers to simulate well the snow 10 thermal profile and the snowpack stratigraphy (Armstrong and Brun, 2008;Vionnet et al., 2012). This configuration allows a good computation of the diurnal cycle through the use of fine top layers, while bottom layers are also sufficiently thin to ensure a good computation of the heat conduction at the snow/soil interface. However, these models were rarely used in global atmospheric, climate, and/or hydrological models due to 15 their high computational costs partly due to the use of a large number of layers. For this reason, the multi-layer snow scheme in ISBA was developed using only three layers representing a good compromise between a reasonable simulation of the snow thermal profile (Boone and Etchevers, 2001) and a low computing time. Today, such computational limitations are less of a constraint and a larger number of layers can 20 be used in this scheme. The number of snow layers in ISBA was increased to 12 with two fine layers at the top and the bottom of the snowpack using the following simple where the constants are defined as: δ 1 = 0.01, δ 2 = 0.05, δ 3 = 0.15, δ 4 = 0.5, δ 5 = 1, δ 9 = 1, δ 10 = 0.5, δ 11 = 0.1, and δ 12 = 0.02 m. For a snow depth below 0.1 m, each layer has the same thickness of 0.00833 m. When the snow depth is above 0.2 m, the thicknesses of the first and the last layers reach their constant values of 0.01 and 0.02 m 10 respectively to reasonably resolve the diurnal cycle and the snow/soil heat exchanges. However, to keep as much as possible the information of an historical snowfall event, the grid thicknesses are updated only if the two first layers or the last layer become too small or too large. This condition can be summed-up as follows: 15 For example, for a total snow depth of 1 m, if the thickness of the top layer becomes lower than 0.005 m or greater than 0.015 m at the beginning of a time step, the layer thicknesses of the entire snowpack are recalculated with Eq. (9) and the snow mass and heat are redistributed accordingly. A similar algorithm was also developed for the 6 and 9 layer cases, but these results are not reported here. In terms of snowpack 20 layering, the main difference with the Crocus scheme is the fact that the total number of layers is constant, while in Crocus only the maximum number of layers is specified 6745 (typically 20 or 50) and the model dynamically uses a number of layers which varies in time within this pre-defined constraint .

Snow compaction and snowdrift events
In the new version of the snow scheme, the evolution of snow density in each layer is due to snow compaction resulting from changes in snow viscosity (Brun et al., 1989) 5 and snow densification of surface layers during snowdrift events (Brun et al., 1997). Snowdrift is assumed to occur when wind velocity exceeds a threshold value that depends on snow surface characteristics. This process is especially important for simulating the evolution of the snow density over polar regions. Brun et al. (1997) pointed out that this process is also critical for reproducing the snow thermal 10 conductivity and the snow temperature profile over these regions. Therefore, the time tendency of snow density in each layer is computed as follows: where ρ drift (kg m −3 ) is the maximum density equal to 350 kg m −3 below which the snow densification occurs during snowdrift events, τ drift (s) the compaction rate of 15 this process (Appendix B), and σ (Pa) the vertical stress in each layer. This stress is computed as the weight of the overlaying layers. At the top of the snow pack, half the mass of the uppermost layer is used. The vertical stress in each layer is then given by: The snow viscosity is a function of snow density, temperature, and liquid water content, W l (kg m −2 ), and it is given as follows: where η 0 (Pa s −1 ) is a reference viscosity equal to 7 622 370 Pa s −1 , ρ 0 (kg m −3 ) is 5 a reference density equal to 250 kg m −3 , W lmax (kg m −2 ) represents the maximum liquid water holding capacity (e.g. Sect. 2.2) and the constants a η = 0.1 K, b η = 0.023 m 3 kg −1 , and ∆T η = 5 K. The viscosity dependence on snow temperature is limited according to Schleef et al. (2014) who pointed out that the impact of snow temperature on snow densification becomes negligible at low temperatures. The last 10 dimensionless function, f W , describes the decrease of viscosity in presence of liquid water. Finally, the snowfall density is computed as previously (Eq. 7).

Transmission of solar radiation and snow albedo
The absorption of incident shortwave solar radiation, R SW (W m −2 ), within the pack is now solved over three spectral bands according to Brun et al. (1992). The first band 15 ([0.3-0.8] µm) represents the ultra-violet and visible range, while the two others ([0.8-1.5] and [1.5-2.8] µm) represent two near-infrared ranges. The total net shortwave radiation, Q sn , absorbed by the snow level i , is the sum of the absorption in each spectral bands, k, and is given by: where ω is the empirical weight of each spectral bands equal to 0. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | coefficient of snow, β sn , depends on density and optical diameter of snow. The snow albedo, α sn , is a function of the snow optical diameter and of the age of the first layer of the snowpack. The age dependency is limited to the first band (visible range) and aims to represent the decrease of the snow albedo by impurities from deposition in a very simple way. Indeed, trace amount of light-absorbing impurities can significantly 5 reduce snow albedo in the visible range but have no effect on the near-infrared range (Warren, 1984). In each band, both the albedo and the extinction coefficient of snow are computed according to Brun et al. (1992) as follows: where A sn is the age of the first snow layer expressed in days, A ref a reference age 15 set to 60 days that modulates the snow albedo decrease due to impurities, P a (Pa) is the near surface atmospheric pressure, and P ref (Pa) a reference pressure equal to 870 hPa. The optical diameter of snow is simply given by Eq. (10) but is now also dependent on snow age: where g 3 is the rate of increase of the optical diameter of snow with snow age. It is set to 0.5 × 10 −4 m day −1 through calibration. The motivation to add this snow age dependency on snow optical diameter is discussed in Sect. 6.

Effects of soil organic carbon on soil hydraulic and thermal properties
North-Eurasian soils are rich in organic carbon as shown in Fig. 1. This figure   5 represents the soil organic carbon content of two soil horizons (0-30 and 30-70 cm) aggregated at a 0.5 • by 0.5 • horizontal resolution and estimated from the Harmonized World Soil Database (HWSD; http://webarchive.iiasa.ac.at/Research/ LUC/External-World-soil-database/HTML/) at a 1 km resolution from the Food and Agricultural Organization (FAO, 2012). The parameterization of the impact of soil 10 organic carbon on hydraulic and thermal properties in ISBA is based on pedotransfer functions of Boelter (1969), and on the work by Letts et al. (2000) and Lawrence and Slater (2008). The pedotransfer functions of Boelter (1969) link the soil water retention at different pressure levels to the fiber content of a peat soil. Letts et al. (2000) describe the vertical profile of hydraulic properties such as soil matric potential and 15 hydraulic conductivity at saturation for a typical organic soil. The hydraulic properties change sharply from the near surface where peat is weakly decomposed (fibric soil) to the sub-surface with moderately and well decomposed peat (hemic and sapric soils respectively). Lawrence and Slater (2008) proposed a linear combination of such soil organic properties with the standard mineral soil properties. 20 In ISBA, before averaging soil organic with mineral properties, a typical peat soil profile is computed for the model soil grid using a power function for each hydraulic property, α peat , found in Table 1. For each soil layer i , this function is described as: where z (m) is the depth of the considered soil grid node, α fibric and α sapric the fibric and sapric parameter values (Table 1), d fibric (m) the depth arbitrarily set to 0.01 m where the profile starts to depart from fibric values, and d sapric (m) the depth of 1 m where the soil properties reach the sapric values according to Letts et al. (2000).
To determine the organic fraction of soil, the density profile of the soil carbon must be 5 known for the entire soil grid. Using the HWSD database, the soil carbon densities in the first 0.3 m, ρ top (kg m −3 ), and the remaining 0.7 m below, ρ sub (kg m −3 ), are known: where S top and S sub (kg m −2 ) are the topsoil and subsoil organic carbon contents respectively, ∆d top and ∆d sub (m) the thicknesses of each observed soil horizon 10 (0.3 and 0.7 m respectively). We extrapolate the density present below 1 m from this observed near-surface profile (Eq. 20). The extrapolation assumes that the carbon profile decreases sharply with soil depth according to a power function. The shape of this function is given by the observed profile if the topsoil organic carbon density is superior to the subsoil density. Otherwise, the density of soil carbon below a 1 m depth, 15 ρ deep (kg m −3 ), is taken equal to the subsoil density: where ∆d deep (m) is an infinite soil thickness taken arbitrarily equal to 1000 m.
Finally, the soil carbon density profile, ρ soc (kg m −3 ), over the entire soil grid is 20 computed using these three soil horizons and a simple linear interpolation at each 6750 grid node that conserves the total soil carbon mass (Fig. 2). The fraction of the soil that is organic, f soc , in each layer is determined assuming this simple relationship: where ρ om (kg m −3 ) is the pure organic matter density equal to 1300 kg m −3 (Farouki, 1986) and w sat, peat the porosity of the peat soil profile computed using Eq. (19) and 5 Table 1. As in Lawrence and Slater (2008), this fraction is used to combine the standard mineral soil properties with soil organic properties using weighted arithmetic or geometric averages, depending on the parameter (Table 1). An example of this method is shown in Fig. 2  The Col de Porte dataset includes many observations at a daily time step for evaluating land surface models. In this study, the observed snow depth, surface albedo and soil temperature at 10 cm are used to evaluate model simulations over the entire period. The snow water equivalent (SWE) is also used for this model evaluation but daily values are only available from 2001 to 2011. Snow depth is measured using ultra-5 sound depth gauges with an accuracy of 1 cm. Surface albedo is computed as the total daily reflected solar flux divided by the total daily incoming solar flux. We estimate the uncertainty in surface albedo to be about 10 % based on the 10 % uncertainty in observed radiative fluxes reported by . Soil temperature is measured using automatic probes with an accuracy of 0.1 K. SWE is measured using a cosmic ray 10 sensor placed on the ground and exhibits an uncertainty of 10 %.Three skill scores are used to compare model results to the observations. The mean annual bias measures the capability of the model to represents the observed mean. To evaluate the model ability to represent the observed day to day variability, two statistical quantities are used; the square correlation (r 2 ), and the centered root mean square error (c-rmse).

15
It is computed by subtracting the simulated and observed annual means from their respective time series before computing a standard root mean square error.

Model configuration
Four simulations were done to evaluate the effect of the different changes in the snow parameterization detailed in Sect. 3:
-SNL is similar to CTL in terms of snow compaction and albedo but uses the new snow layering with 12 snow layers described in Sect. 3.1.1. -CPT uses 12 snow layers as in SNL but the compaction and the snow densification of surface layers during snowdrift events are computed using formulations of Brun et al. (1989Brun et al. ( , 1997, both described in Sect. 3.1.2. -NEW uses all the package of snow equations described in Sect. 3.1: 12 snow layers, the new snow compaction/densification, but also the spectral 5 representation of the snow albedo (Sect. 3.1.3).
For all of the simulations, the snow is assumed to cover the entire grid-cell (the snow fraction set to 1) as long as the snow remains present. The effective roughness length of snow is set to its usual value of 0.001 m. The grid-cell is assumed to be entirely covered by grass with a root depth of 1 m, the leaf area index varies from 0.1 10 in winter to 1 in summer, and the snow-free surface albedo is prescribed as 0.2. The model calculates soil temperature, moisture and ice content in each of the 14 soil layers corresponding to a soil depth of 12 m. The model was run with a 15 min time step from 1 August 1993 to 31 July 2011. The model was spun-up by performing fifty iterations of the first two years (August 1993 to July 1995). This spin-up represents a total of one 15 hundred years, and this was determined to guarantee that the water and heat profiles were equilibrated over the 12 m soil depth of ISBA. Results are then evaluated over the entire period. allow a discrimination between the four simulations, even if the c-rmse of the NEW experiment is the best for seven of the ten years. The surface albedo from the NEW simulation is clearly better than the albedo from the other experiments: bias and crmse are the best for all years (Fig. 4). The soil temperature bias and c-rmse are also reduced by the NEW experiment (for ten of seventeen years) compared to the 5 other simulations. Thus, accounting for different spectral bands within the snow albedo calculation has a significant positive impact on the energy balance of the snow-soil system. The average seasonal cycle of snow depth, SWE, surface albedo and soil temperature at 10 cm represented in Fig. 5 highlights the qualities and weaknesses 10 of the different parameterizations by focusing on the snow season (October-May). The corresponding statistics for the winter (DJF), spring (MAM) and the entire period are given in Table 2. The comparison of SNL-CTL indicates that the increase in number of snow layers from 3 to12 improves the snow depth, SWE and winter soil temperature simulation. Change in snow compaction (SNL-CPT) improves the seasonal cycle 15 of snow depth and SWE and especially the maximum value. The seasonal and global biases in Table 2 verify this result and show the same behavior for winter soil temperature, although it is difficult to see visually from Fig. 5. For these three variables, the simulated time variability is also improved from CTL to SNL to CPT as shown by the other seasonal and global scores (c-rmse and r 2 ) in Table 2. Finally, the new 20 spectral albedo scheme (CPT-NEW) has a drastic impact on the snowpack simulation in spring. As shown by Fig. 5 and Table 2, the new spectral albedos clearly improve the simulation of other variables during this period. They induce a sharp springtime snowmelt with a strong decrease in snow depth and SWE. The snow insulation during spring is thus less important and allows the soil surface to warm up faster. As a result, 25 the model is capable of reproducing the strong soil warming observed in April (Fig. 5).

Figures 3 and 4 show an overview of the four simulations performed at the Col
Not surprisingly, the soil temperature skill scores for spring and the whole period are drastically improved although there is a slight degradation in winter.  Figure 6 shows daily mean time series of the snow density and temperature profiles averaged over the snow season for each experiment. With only 3 snow layers (CTL), the density distribution is more uniform than using the new snow layering scheme with 12 layers (SNL). The significant densification of the bottom layers in SNL is the main process responsible for the snow depth and SWE improvements observed in Fig. 5   5 and Table 2. In terms of snow temperature, SNL induces warmer bottom snow layers in winter due to the densification of the lowest snow model layers and the insulation from other overlying layers. This explains the skill scores improvement found in winter soil temperature in Table 2. The new snow compaction scheme (CPT) tends to increase the density contrast between the top and the bottom snow layers. The snowpack is also denser than with SNL leading to the strong decrease in snow depth observed in Fig. 5 and to the better skill scores in snow depth and SWE over each period (Table 2).
CPT also results in a small warming at the bottom of the snowpack which slightly heats the soil temperature compared to SNL. Finally, the spectral albedo scheme (NEW) has a limited effect on the snow density profile but results in a slightly colder 15 snowpack than in CPT and even SNL (not show) due to the large daily winter albedos seen in Fig. 5. This is the main reason for the lower winter soil temperatures with NEW than CPT and SNL (Table 2). 20 The experimental design used here is close to that proposed by Brun et al. (2013) Dee et al. (2011) and an evaluation of its performance is provided in Berrisford et al. (2011). For precipitation, the monthly ERA-I precipitation are rescaled to match the observed Global Precipitation Climatology Center (GPCC) Full Data Product V5 (http://gpcc.dwd.de) as proposed by Decharme and Douville (2006a). This method 5 conserves the 3-hourly chronology of the ERA-I precipitation but ensures a reasonable monthly amount (Szczypta et al., 2012). Brun et al. (2013) pointed out the significantly better performance of this ERA-I scaled GPCC forcing product in simulating North-Eurasian snowpack variables compared to the ERA-I precipitation or other "state of the art" global scale atmospheric forcings. 10 To evaluate snow and soil temperature simulations, several in situ dataset are used. As in Brun et al. (2013), the Historical Soviet Daily Snow Depth (HSDSD; http: //nsidc.org/data/docs/noaa/g01092_hsdsd/index.html) compiled by Amstrong (2001) was used in the current study. It consists in daily snow depth measurements taken at synoptic stations following the World Meteorological Organization (WMO) standards. 15 WMO requires the measurements to be taken in bare ground open areas or clearings with regular grass cutting. These snow depth data are therefore representative of open areas of bare ground or those covered with very short grass. This dataset starts in 1881 with a few stations and ends in 1995. Considering that ERA-I starts in 1979, the model simulations are done from 1979 to 1993 according to Brun et al. (2013). 263 20 HSDSD stations are available over this period with approximately half of them without any missing data. We chose to use only the stations where the difference between the local and the ERA-I elevation is less than 100 m to avoid temperature biases for instance that would be directly due to the low resolution of ERA-I. We also only kept the stations where the number of days with a non zero snow depth measurement over 25 the entire period is superior to 100 days, and that have at least 8 days with snow measurement per year. With this filter, the number of available stations decreases to 158, which remains acceptable. Most stations are located in Russia and Western-Siberia with only a few in Eastern-Siberia (Fig. 7). The second source of observations is the Russian Historical Soil Temperature (RHST) dataset compiled by Zhang et al. (2001) over Siberia (http://data.eol.ucar. edu/codiac/dss/id=106.ARCSS078). Data coverage extends from the 1800s through 1990, but is not continuous. We compared our model results over the 1979-1990 period. As for snow depth, soil temperature stations are subject to WMO standards 5 and are located in open area sites. We used the same criteria as for snow depth. Only stations with local elevations close to the ERA-I altitude (less than 100 m difference) are used. In addition, only stations with at least 36 months of observations (at least 3 years out of 12) are kept. Most soil temperature sites are collocated with snow depth sites (Fig. 7). Measurements were taken at 20, 80, 160 and 320 cm depth. For each 10 depth, 95, 48, 48, and 82 stations, respectively, were available for model evaluation. The spatial distribution of these stations is shown in Fig. 7 for 20 and 160 cm depths.

Numerical experiment design and observational dataset
To quantify the capability of the model to simulate the permafrost characteristics, two datasets are used. The first dataset is the Circum-Artic Map of Permafrost and Ground Ice Conditions (http://nsidc.org/data/ggd318) edited by Brown et al. (2002). This 15 dataset is available at a 0.5 • ×0.5 • resolution and shows the continuous, discontinuous, isolated and sporadic permafrost boundaries. The second dataset is an estimate of the active layer thickness over North-West-Siberia. This data set is based on the map of landscapes and permafrost conditions in Yakutia (http://doi.pangaea.de/10.1594/ PANGAEA.808240). It gives access to the mean and standard deviation of the most 20 probable active layer thickness in each grid box at 0.5 • × 0.5 • resolution. All details can be found Beer et al. (2013).

Model configuration
Three experiments using the ISBA land surface model forced by the ERA-I scaled GPCC atmospheric dataset are performed using the same configuration. In addition 25 to the CTL (old snow scheme) and NEW (new snow scheme) experiments already described in Sect. 4, we performed one simulation using the parameterization of the impact of the soil organic carbon on the hydrologic and thermal soil properties. This 6757 last experiment, called NEW-SOC, uses the new snow and soil-property schemes described in Sects. 3.1 and 3.2, respectively. As previously, the model determines the temperature, liquid water and ice content evolution in each of the 14 soil layers corresponding to a total soil depth of 12 m. The model is run with a 15 min time step from 1 January 1979 to 31 December 1993. The model's spin-up uses twenty iterations of the first five years (1979)(1980)(1981)(1982)(1983) of the atmospheric forcing, representing a total of one hundred years. In ISBA, we use a series of twelve sub-grid independent patches per grid cell in order to account for land cover heterogeneity. Land cover parameters such as leaf area index (LAI), vegetation height, vegetation/soil albedos, and rooting depth are prescribed for each sub-grid patch. The dominant patches present in the model over the Northern-Eurasian region are bare soil, grassland/tundra, deciduous forest, coniferous boreal forest, and C3 crops in the South. The fraction of each surface type within each grid box is used to compute the grid box average of the water and energy budgets. Some other processes, such as surface runoff, dripping from the canopy reservoir, and 15 soil infiltration account for sub-grid parameterizations. More details can be found in Decharme and Douville (2006b) and Decharme et al. (2013).
For all of the simulations, the grid-cell fraction covered by snow evolves according to the simulated snow depth and is different for bare soil and vegetated areas (Appendix C) in each land cover patch. As was the case for the Col de Porte 20 experiment, the effective roughness length of snow retains its usual value of 0.001 m. The land surface parameters used by ISBA are specified according to the 1 km resolution ECOCLIMAP-II database (Faroux et al., 2013). LAI, vegetation height, and vegetation/soil albedos are prescribed for the twelve vegetation sub-grid patches based on a mean annual cycle at a 10 day time step. The rooting depth is specified for 25 each vegetation type according to Canadell et al. (1996). It ranges from 0.5 to 1.5 m for tundra and temperate grassland, and from 2 to 3 m for forest. The soil textural properties are given by the HSWD database at 1 km resolution. The topographic information is specified according to the 30-arcsecond resolution GTOPO30 data set, and the topographic indexes used in the TOPMODEL runoff parameterization are given by the 1 km resolution HYDRO1K product (http://eros.usgs.gov/#/Find_Data/Products_ and_Data_Available/gtopo30/hydro). Figure 7 presents a first quantitative comparison between the observed and simulated 5 snow depth and soil temperature over Northern-Eurasia. Because in situ observations were collected in bare ground open areas and/or clearings with regular grass cutting following the WMO standards as mentioned previously, they are compared to snow depths and soil temperature profiles simulated by the ISBA bare soil sub-grid patch alone. This patch exhibits conditions which are closest to those at the corresponding 10 field sites, as is generally the case for ISBA in this kind of comparison . The simulation represented here is the NEW-SOC experiment that seems to capture well the snow depth and soil temperature spatial distributions. For snow depth, the latitudinal gradient is well respected. The lower soil temperature along a southwest-northeast transect is also well simulated. 15 The seasonal cycles of daily snow depths and monthly soil temperatures (Fig. 8) clearly show the biases of the CTL simulation and the improvements due to the new snow and soil representations. The seasonal cycles and the global skill scores are computed using the measurements and simulations for all stations over the entire observed periods. ISBA globally underestimates the snow depth from December 20 through February with no clear difference between CTL and NEW (or NEW-SOC). However, the springtime snow ablationis drastically improved by the new snow scheme inducing a better simulated seasonality. This fact is confirmed by some other quantitative comparisons. The average number of days per year with observed snow on the ground for all in situ stations is 150.7 days. CTL simulates 158.7 days against 25 151.5 days for NEW. On average, the last day of the snow season is day number 281.6 when starting on July first. CTL goes beyond this date by more than 9 days while for NEW it is only 2 days (day number 283). Theses results are consistent withthe model 6759 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | evaluation at the Col de Porte field site (Sect. 4). As could be expected also, the new physical soil properties (NEW-SOC) play a minimal role in the snow depth simulation. The seasonal cycle of the soil temperature profile confirms that the new snow scheme induces a warmer soil in winter compared to CTL, and it strongly reduces the cold bias of CTL. The effect of soil organic carbon is especially observable during spring and 5 summer. NEW exhibits a warm bias for each soil horizon while NEW-SOC, with more insulating soils, reduces this weakness. These improvement in snow depth and soil temperature are confirmed by the spatial distributions of their seasonal skill scores (bias and c-rmse). Figure 9 shows the spatial distributions of snow depth seasonal skill scores (bias and c-rmse) during winter and 10 spring. No clear discrepancies among these simulations appear in winter while the bias and c-rmse of many stations are improved in spring by the new snow scheme. The springtime snow depth is simulated in an acceptable manner by NEW, while CTL exhibits a significant overestimation. This fact is confirmed by global scores given in each of the panels. In winter, regardless of the experiments, ISBA underestimates 15 snow depth measurements at many stations, especially in the Northern and Western parts of the domains. The spatial distribution of soil temperature seasonal skill scores simulated at 20 and 160 cm depth during winter is given in Fig. 10. Regardless of the region, the generalized cold bias found over all stations with CTL is drastically reduced with the new snow scheme and the interannual variability (c-rmse) is largely 20 improved. In summer, as was already shown in Fig. 8, NEW-SOC is in better agreement with observations regardless of the soil horizon even if a slight cold bias appears at the subsurface. (Fig. 11). The NEW experiment overestimates the temperature profile measurements at many stations near the surface, but less-so at a 320 cm depth. So, it seems that the subsurface cooling in the NEW-SOC experiment is too intensive. But in Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | are in fact due to error compensation between the cold and warm biases simulated in the western and eastern part of the domain, respectively. The effect of soil organic carbon content on soil temperature profile is also especially observable in terms of the simulated permafrost characteristics. The observed and simulated locations of permafrost boundaries are compared in Fig. 12. Regardless of 5 the experiment, ISBA generally simulates acceptable boundaries even if the permafrost limit extends slightly too far south in the western part of the domain. This figure also shows the spatial distribution of active layer thicknesses simulated by the NEW and the NEW-SOC experiments. The active layer thickness in the model is computed as the maximum depth reached each year by the 0 • C isotherm in the soil approximated 10 via a linear extrapolation between the last node from the surface where the soil temperature is positive and the first where it is negative. As expected from the lower summer soil temperatures with NEW-SOC (Figs. 9 and 11), the active layer is shallower. However, this comparison with the limits of different permafrost types does not allow the determination of which simulation leads to the most accurate active layer 15 thicknesses, unless it is considered that the continuous and discontinuous permafrost exhibit generally active layers shallower than 1.5 or 2 m. With this very subjective criterion, NEW-SOC seems to simulate a more accurate spatial distribution of the active layer thickness. To partially verify this hypothesis, Fig. 13 shows the estimated and simulated active layer thicknesses over the Yakutia region. Estimations from Beer 6 Discussion and conclusion 25 In this study, the impact of improved representation of snowpack and soil properties in the ISBA LSM to simulate snow characteristics and soil temperature profiles over cold regions was analysed. ISBA's representations of snow layering, albedo, and compaction were updated by incorporating some parameterizations of the detailed Crocus snowpack model. In addition, a simple parameterization of the soil organic carbon effect on hydraulic and thermal soil properties was introduced based on previous work (Boelter, 1969;Letts et al., 2000;Lawrence and Slater, 2008). The model is evaluated first over the Col de Porte field site in the French Alps  in order to isolate the changes in the snowpack parameterization, and second over the North-Eurasian region to analyze the model's ability to simulate snow depth, soil temperature profile and permafrost characteristics. Changes in the snowpack parameterizations induce noticeable improvements in the 10 simulated snow depth, SWE, surface albedo and soil temperature at the Col de Porte (field) site. The new snow layering algorithm with 12 layers permits a refinement of the vertical distribution of density and temperature in the snowpack leading to slight improvements in simulated snow depth, SWE, and soil temperature during winter. The densification of the snowpack with the new compaction scheme, which increases the 15 density contrast between the top and the bottom snow layers, has a significant positive impact on snow depth and winter soil temperature. Finally, the new spectral albedo scheme clearly improves the simulation of the springtime surface albedo that allows a better simulation of the snowpack characteristics and soil temperature during ablation at the end of the snow season. 20 It must be noted that the large improvement in snow albedo in spring is mainly due to the use of snow age in the diagnostics of the optical diameter of snow (Eq. 16). Without this parameterization, the surface albedo is strongly overestimated in winter and, to a lesser extent in spring at the Col de Porte field site, with a larger bias and c-rmse for all variables compared to the new version of ISBA (not shown). The optical 25 diameter of snow strongly controls the near-infrared albedo, while impurities mostly affect the albedo in the visible spectrum (Wiscombe and Warren, 1981). This increase of snow optical diameter with time is necessary to represent well the decrease in spectrally integrated albedo with age. However, the increase of snow optical diameter is not only a function of snow density as parameterized by Anderson (1976) in Eq. (8), but it is also due to snow metamorphism, which is macroscopically driven by snow temperature and snow thermal gradients. Several complex parameterizations exist to explicitly represent the evolution of snow optical diameter according to these processes (e.g. Carmagnola et al., 2014). Nevertheless, for the sake of simplicity, we just use 5 a snow age dependency in the diagnostic of snow optical diameter with a limitation at fifteen days (Eq. 16). This simple diagnostic allows the model to reasonably match the explicit computation of the optical diameter of snow simulated in the Crocus model (not shown). The good results of the ISBA model at the Col de Porte field site reinforce this choice. 10 The positive impacts of the new ISBA snow scheme are confirmed when tested over the North-Eurasian region with an important number of open field in situ snow depth and soil temperature stations. Winter snow depths are slightly better simulated with the new version and the winter soil temperature cold bias obtained with the old version of ISBA is clearly reduced. This fact confirms that the physics used in snow schemes is of 15 primary importance for adequately simulating the snow insulating effect that prevents soil from getting too cold in winter (Slater et al., 2001;Luo et al., 2003;Gouttevin et al., 2012;Paquin and Sushama, 2015). Another important impact of changes in the ISBA snow scheme over the North-Eurasian region is seen in spring when the snowmelt is well reproduced. As shown over the Col de Porte (field) site, this is mainly 20 due to the new parameterization of spectral snow albedo.
Nevertheless, regardless of the model version used, simulated winter snow depths are generally underestimated compared to in situ observations. The cause of this underestimation is not trivial. The first source of uncertainty can be attributable to the GPCC precipitation measurements that do not account for wind undercatch leading to 25 a possible underestimation of solid precipitation during winter (Adam and Lettenmaier, 2003;Brun et al., 2013). Besides uncertainties related to the atmospheric forcing, the snow depth underestimation can be due to the non-explicit representation of snow metamorphism. Indeed, in similar experimental conditions over the Northern Eurasian region, the winter snow depth simulated by the detailed Crocus snowpack model did not exhibit the same problem (Brun et al., 2013) and the main remaining difference between Crocus and ISBA is now restricted almost entirely to the explicit simulation of snow metamorphism. In Crocus, the viscosity of layers composed of facetted crystals and depth hoar snow types is increased , which leads to reducing 5 the overall compaction rate of snowpack undergoing temperature conditions conducive to such snow types, and this is consistent with the situation described above. Taking into account soil organic carbon in soil physical properties logically plays a minimal role in the simulated snowpack behaviour. However, this process has drastic impacts on the summer soil temperature profile because it allows the soil to remain cool during spring and summer as shown in previous studies (Bonan and Shugart, 1989;Lawrence and Slater, 2008;Dankers et al., 2011). Consequently, the spatial distribution of the permafrost active layer thickness simulated by the new version of ISBA is in better agreement with estimations from Beer et al. (2013) over the Yakutia region. This result is in agreement with Paquin and Sushama (2015) who showed that the hydraulic 15 and thermal properties of soil organic carbon partly control the thickness of the active layer during summer. However, spatial observations of permafrost characteristics on the global scale are still very scarce, and if available, they are static and don't allow the study of long term trends and inter-annual variability.
This model validation should ideally be extended over all cold regions (e.g. North

20
America, Greenland, etc.) but considering that North-Eurasia is representative of such regions, some important conclusions are confirmed by this study: -An adequate simulation of snow layering and snow compaction/densification is important in order to represent well winter snowpack characteristics and the soil temperature profile. -To account for soil organic carbon in terms of the soil physical properties drastically impacts the simulation of summer soil temperture profile and hence the permafrost active layer tichkness and its spatial distribution.
Finally, these conclusions underscore the fact that the representation of snowpack characteristics and soil thermal processes are of primary importance for studying 5 permafrost vulnerability under climate change conditions, especially if the continental carbon cycle is considered due to the strong interaction between soil carbon degradation and soil thermal processes.

Appendix A: A snow thermal conductivity
The snow thermal conductivity is computed as a function of snow density following Yen 10 (1981). It also accounts for vapor transfer in the snow using a simple parameterization from Sun et al. (1999). This process is especially important at low snow densities and at high altitude. So the snow thermal conductivity, λ sn (W m −1 K −1 ), in each layer is given by: where λ ice (W m −1 K −1 ) is the thermal conductivity of ice equal to 2.2 W m −1 K −1 , ρ w (kg m −3 ) the water density, P a (Pa) the air pressure, P 0 (Pa) a reference pressure equal to 1000 hPa, and the coefficients k 1 = −0.06023 W m −1 K −1 , k 2 = 2.5425 W m −1 and k 3 = 289.99 K. the potential for snow erosion for each snow layer is computed as a function of snow density: where ρ sn min = 50 kg m −3 is the minimum density of snow, ρ mob a reference density of 295 kg m −3 , and the dimensionless constant a mob = 1.25. Secondly, a driftability index, 5 Γ drift , combining the mobility index and the near surface atmospheric wind speed: where κ v = 1.25 is a dimensionless coefficient for gust diagnosis from average wind speed, and the constants a Γ = 2.868 and b Γ = 0.085 s m −1 . A positive value of Γ drift indicates that snowdrifting can occur. Compaction rate from the surface is then 10 propagated to the layers beneath, following an exponential decrease, until it meets a snow layer having a negative driftability index. For each layer, this compaction rate is computed as follows: where π τ (s) is a time constant of one day, and the constants a τ = 10 and b τ = 3. 25 0.8 m for grassland/tundra. Finally, the height of crop types is related to an exponential function of LAI and has a height of 1 m before maturity defined as a LAI of 3.5 m 2 m −2 . More details on these physiographic parameters can be found in Masson et al. (2003).   Data Discuss., 5, 29-45, doi:10.5194/essdd-5-29-2012 Figure 2. Parameterization of the effect of soil organic carbon (SOC) on soil hydraulic and thermal properties. The soil organic carbon density profile, ρ soc is given by Eq. (21) using a top soil organic carbon content of 10 kg m −2 , a sub soil content of 15 kg m −2 , and via a simple linear interpolation at each soil grid nodes that conserves the total soil carbon mass. The fraction of the soil that is organic, f soc , in each layer is determined assuming a simple relationship between this last soil organic carbon density profile and an idealized peat soil density profile (Eq. 22). Examples for the soil porosity, w sat , the soil saturated hydraulic conductivity, k sat , and the soil heat capacity, c, are given. Dotted lines represent vertical homogeneous mineral soil properties, dashed lines the idealized peat soil properties, and plain lines the resulting combined soil properties using averaging method sums-up in Table 1 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Figure 5. Daily mean annual cycles of snow depth, SWE, surface albedo, and soil temperature at 10 cm depth simulated (colours) and observed (black) at the Col de Porte field site. The corresponding skill scores are given in Table 2. Over all panels, the grey shadow corresponds to the uncertainty in in situ measurements as discussed in Sect. 4.1. The observed snow depth exhibits an accuracy of ±1 cm, the soil temperature is measured with a precision of ±1 K, while uncertainties in SWE and surface albedo is near ±10 %.