Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model

In this study we analyzed how an improved representation of snowpack processes and soil properties in the multilayer snow and soil schemes of the Interaction SoilBiosphere-Atmosphere (ISBA) land surface model impacts the simulation of soil temperature profiles over northern Eurasian regions. For this purpose, we refine ISBA’s snow layering algorithm and propose a parameterization of snow albedo and snow compaction/densification adapted from the detailed Crocus snowpack model. We also include a dependency on soil organic carbon content for ISBA’s hydraulic and thermal soil properties. First, changes in the snowpack parameterization are evaluated against snow depth, snow water equivalent, surface albedo, and soil temperature at a 10 cm depth observed at the Col de Porte field site in the French Alps. Next, the new model version including all of the changes is used over northern Eurasia to evaluate the model’s ability to simulate the snow depth, the soil temperature profile, and the permafrost characteristics. The results confirm that an adequate simulation of snow layering and snow compaction/densification significantly impacts the snowpack characteristics and the soil temperature profile during winter, while the impact of the more accurate snow albedo computation is dominant during the spring. In summer, the accounting for the effect of soil organic carbon on hydraulic and thermal soil properties improves the simulation of the soil temperature profile. Finally, the results confirm that this last process strongly influences the simulation of the permafrost active layer thickness and its spatial distribution.


19
Snowpack properties are known to be of primary importance for understanding the water 20 and energy budgets of the land surface, especially in mountainous and boreal regions. From 21 autumn to spring, solid precipitation is stored within the snowpack thereby modifying the 22 terrestrial albedo and roughness length, and impacting the radiative and energy fluxes at the 23 soil/atmosphere interface. During spring, the fresh water released by snowmelt contributes to soil 24 infiltration, intense streamflow and large seasonal flood events, and it directly modulates the land 25 surface evapotranspiration [Poutou et al. 2004;Niu and Yang 2006;Decharme and Douville 26 2007]. Snowpack also acts as an insulating layer at the surface which prevents significant heat loss 27 in the winter. Over North-Eurasian regions, as discussed by Paquin and Sushama [2015], this last 28 process controls the temperature of the permafrost. It is defined as a soil that remains below 0°C 29 for two or more consecutive years, and it has a significant influence on the summer permafrost 30 active layer thickness, defined as the maximum annual thaw depth. In summary, snowpack 31 properties drastically influence soil/atmosphere interactions during a large part of the year through 32 their impacts on many land surface processes. 33 Beside the importance of snowpack properties for understanding the water and energy 34 budgets of the land surface in northern regions, the physical properties of soil organic carbon (or 35 peat soil) also play a significant role. North-Eurasian soils are very rich in organic carbon because 36 the low soil temperatures in this region inhibit decomposition of dead plant material that 37 accumulates over time, thereby forming peat deposits. Soil organic carbon exhibits very different 38 hydraulic and thermal properties than mineral soil [Boelter 1969;Letts et al. 2000]. It is 39 characterized by a very high porosity, a weak hydraulic suction, and a sharp vertical hydraulic 40 conductivity profile from high values at the surface to very low values at the subsurface. This 41 generally induces a relatively wet soil with a shallow water table [Letts et al. 2000]. Its low 42 thermal conductivity and its relatively high heat capacity act as an insulator for soil temperature 43 that prevents the soil from significant warming during the summer [Bonan and Shugart 1989; where the soil matric potential ψ is defined using Equation 1. Thus, this scheme induces 161 dependencies of water phase changes to soil textures and to the degree of soil humidity. The 162 coarser the soil texture, the larger the quantity of water that will freeze at a given temperature. As 163 the soil becomes dry, the temperature that allows freezing drops. More

Snowpack internal processes 167
The original ISBA explicit multi-layer snow scheme developed by Boone and Etchevers 168 [2001] is a snowpack scheme of intermediate complexity made in order to take into account for 169 some processes such as snow mass and heat vertical redistribution, snow compaction, water 170 percolation and refreezing, and explicit heat conduction at the snow/soil interface. Many of theses 171 processes, such as snow compaction or absorption of solar energy, are based on works of 172 Anderson [1976] and Loth et al. [1993]. The thermal conductivity of snow (Appendix A) is 173 computed via the snow density [Yen 1981]. An additional term depends on the snow temperature 174 to account for vapor transfer through the snowpack [Sun et al. 1999]. The time evolution of the 175 snow mass is linked to snowmelt, water freezing, evaporation, and liquid flow. The liquid water 176 content into the snowpack is simulated as a succession of bucket-type reservoirs. A maximum 177 liquid water holding capacity (W lmax ) is computed in each layer. It varies from 3% to 10% of the 178 snow mass according to a decrease in snow density after Anderson [1976]. A liquid water flux is 179 generated when the liquid water content exceeds this threshold. More details can be found in 180 Boone and Etchevers [2001] and only internal physical processes of the snowpack discussed in 181 this study are described below. 182

Snow layering 183
In the original ISBA explicit snow scheme, three-layers are used for snow layering because 184 it is considered to be the minimum number required to resolve adequately the snow thermal profile 185 within the snowpack [Lynch-Stieglitz 1994; Loth and Graf, 1998; Boone and Etchevers 2001]. The 186 algorithm that computes the snow grid thicknesses ∆z of each layer, i, is described as follows: 187 where h sn (m) is the total snow depth. As long as the snow remains below 0.2m, the fraction of the 189 total depth that defines the thickness of each layer remains with a fine resolution at the top and the 190 base of the snowpack. When the snow depth exceeds 0.2m, the thickness of the first layer remains 191 equal to 0.05m, in order to adequately solve the diurnal cycle of the surface energy balance. In 192 addition, for large snow depth values, the second layer thickness cannot exceed 0.5m because 193 density and heat vertical gradients are generally the largest near the top of the snowpack. The 194 vertical grid is updated at the beginning of each time step before the computation of the other 195 snowpack internal processes. 196

Transmission of solar radiation and Snow albedo 212
The absorption of incident shortwave solar radiation, R SW (W.m -2 ), within the snowpack is 213 solved over a single spectral band. It uses an exponential decrease of incoming radiation with 214 snow depth [Anderson 1976;Loth et al. 1993]. So, the net shortwave radiation Q sn (W.m -2 ) 215 absorbed by the snow level, i, is given by: 216 where α sn is the dimensionless snow albedo, and β sn (m -1 ) the extinction coefficient of snow which 218 is given by : 219 As shown by Bohren and Barkstrom [1974], this extinction of snow is directly related to its 221 density, the optical diameter d opt (m), and a constant C ν = 3.8 10 -3 m 5/2 .kg -1 . The optical diameter is 222 empirically linked to the snow density following a simple polynomial regression established by 223 Anderson [1976]: 224 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, 226 and g 2 = 1.1 10 -13 m 13 .kg -4 were calibrated by Anderson [1976]. The time evolution of snow albedo 227 is modelled in a simple way using time constants after Douville et al. [1995]. A linear decrease 228 rate is used for dry snow and an exponential decrease is used for wet snow while the snow albedo 229 increases linearly with snowfall intensity [Boone and Etchevers 2001]. The snow albedo is 230 constrained to be between its minimum value, α min = 0.5, and its maximum, α max = 0.85. 231 while bottom layers are also sufficiently thin to ensure a good computation of the heat conduction 238 at the snow/soil interface. However, these models were rarely used in global atmospheric, climate, 239 and/or hydrological models due to their high computational costs partly due to the use of a large 240 number of layers. For this reason, the multi-layer snow scheme in ISBA was developed using only 241 three layers representing a good compromise between a reasonable simulation of the snow thermal 242 profile [Boone and Etchevers 2001] and a low computing time. Today, such computational 243 limitations are less of a constraint and a larger number of layers can be used in this scheme. The 244 number of snow layers in ISBA was increased to 12 with two fine layers at the top and the bottom 245 of the snowpack using the following simple algorithm: 246

Changes in explicit snow and soil schemes
where the constants are defined as: δ 1 = 0.01m, δ 2 = 0.05m, δ 3 = 0.15m, δ 4 = 0.5m, δ 5 = 1m, δ 9 = 248 1m, δ 10 = 0.5m, δ 11 = 0.1m, and δ 12 = 0.02m. For a snow depth below 0.1m, each layer has the 249 same thickness of 0.00833m. When the snow depth is above 0.2m, the thicknesses of the first and 250 the last layers reach their constant values of 0.01m and 0.02m respectively to reasonably resolve 251 the diurnal cycle and the snow/soil heat exchanges. However, to keep as much as possible the 252 information of an historical snowfall event, the grid thicknesses are updated only if the two first 253 layers or the last layer become too small or too large. This condition can be summed-up as 254 follows: 255 For example, for a total snow depth of 1m, if the thickness of the top layer becomes lower than 257 0.005m or greater than 0.015m at the beginning of a time step, the layer thicknesses of the entire 258 snowpack are recalculated with Equation 11 and the snow mass and heat are redistributed 259 accordingly. A similar algorithm was also developed for the 6 and 9 layer cases, but these results 260 are not reported here. In terms of snowpack layering, the main difference with the Crocus scheme 261 is the fact that the total number of layers is constant, while in Crocus only the maximum number of 262 layers is specified (typically 20 or 50) and the model dynamically uses a number of layers which 263 varies in time within this pre-defined constraint [Vionnet et al 2012]. 264

Snow compaction 265
In the new version of the snow scheme, the evolution of snow density in each layer is due 266 to snow compaction resulting from changes in snow viscosity [Brun et al 1989] and wind-induced 267 densification of near surface snow layers [Brun et al. 1997]. This wind-driven compaction process 268 is assumed to occur when wind velocity exceeds a threshold value that depends on snow surface 269 characteristics. This process is especially important for simulating the evolution of the snow 270 density over polar regions. Brun et al. [1997] pointed out that this process is also critical for 271 reproducing the snow thermal conductivity and the snow temperature profile over these regions. 272 Therefore, the time tendency of snow density in each layer is computed as follows: 273 where ρ wmax (kg.m -3 ) is the maximum density equal to 350 kg.m -3 below which the snow 275 densification occurs during wind-driven compaction, τ w (s) the compaction rate of this process 276 (Appendix B), and σ (Pa) the vertical stress in each layer. This stress is computed as the weight of 277 the overlaying layers. At the top of the snow pack, half the mass of the uppermost layer is used. 278 The vertical stress in each layer is then given by: 279 The snow viscosity is a function of snow density, temperature, and liquid water content, W l 281 (kg.m -2 ), and it is given as follows: 282 where η 0 (Pa s) is a reference viscosity equal to 7622370 Pa s, ρ 0 (kg.m -3 ) is a reference density 284 equal to 250kg.m -3 , W lmax (kg.m -2 ) represents the maximum liquid water holding capacity (e.g. 285 section 2.2) and the constants a η =0.1K -1 , b η = 0.023 m 3 .kg -1 , and ∆T η = 5K. The viscosity 286 dependence on snow temperature is limited according to Schleef et al. [2014] who pointed out that 287 the impact of snow temperature on snow densification becomes negligible at low temperatures. 288 The last dimensionless function, f W , describes the decrease of viscosity in presence of liquid water. 289 Finally, the snowfall density is computed as previously (Equation 7). 290

Transmission of solar radiation and Snow albedo 291
The absorption of incident shortwave solar radiation, R SW (W.m -2 ), within the pack is now 292 solved over three spectral bands according to Brun et al. [1992]. The first band ([0.3-0.8] µm) 293 represents the ultra-violet and visible range, while the two others ([0.8-1.5] µm and [1.5-2.8] µm) 294 represent two near-infrared ranges. The total net shortwave radiation, Q sn , absorbed by the snow 295 level i, is the sum of the absorption in each spectral bands, k, and is given by: 296 where ω is the empirical weight of each spectral bands equal to 0.71, 0.21, and 0.08 for [0.3-0.8], 298 [0.8-1.5] and [1.5-2.8] µm, respectively. As previously, the extinction coefficient of snow, β sn , 299 depends on density and optical diameter of snow. The snow albedo, α sn , is a function of the snow 300 optical diameter and of the age of the first layer of the snowpack. The age dependency is limited to 301 the first band (visible range) and aims to represent the decrease of the snow albedo by impurities 302 from deposition in a very simple way. Indeed, trace amount of light-absorbing impurities can 303 significantly reduce snow albedo in the visible range but have no effect on the near-infrared range 304 [Warren 1984]. In each band, both the albedo and the extinction coefficient of snow are computed 305 according to Brun et al. [1992] as follows: 306 where A sn is the age of the first snow layer expressed in days, A ref a reference age set to 60 days 308 that modulates the snow albedo decrease due to impurities, P a (Pa) is the near surface atmospheric 309 pressure, and P ref (Pa) a reference pressure equal to 870hPa. The optical diameter of snow is 310 simply given by Equation (10) but is now also dependent on snow age: 311 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 313 m.day -1 through calibration. The motivation to add this snow age dependency on snow optical 314 diameter is discussed in section 6. In ISBA, before averaging soil organic with mineral properties, a typical peat soil profile is 339 computed for the model soil grid using a power function for each hydraulic property, α peat , found 340 in Table 1. For each soil layer i, this function is described as: 341 where z (m) is the depth of the considered soil grid node, α fibric and α sapric the fibric and sapric 343 parameter values (Table 1) The extrapolation assumes that the carbon profile decreases sharply with soil depth according to a 354 power function. The shape of this function is given by the observed profile if the topsoil organic 355 carbon density is superior to the subsoil density. Otherwise, the density of soil carbon below a 1m 356 where ∆d deep (m) is an infinite soil thickness taken arbitrarily equal to 1000m. 359 Finally, the soil carbon density profile, ρ soc (kg.m -3 ), over the entire soil grid is computed 360 using these three soil horizons and a simple linear interpolation at each grid node that conserves 361 the total soil carbon mass (Figure 2). The fraction of the soil that is organic, f soc , in each layer is 362 determined assuming this simple relationship: 363 where ρ om (kg.m -3 ) is the pure organic matter density equal to 1300 kg.m -3 [Farouki 1986] and 365 w sat,peat the porosity of the peat soil profile computed using Equation 19 and Table 1. As in 366 Lawrence and Slater [2008], this fraction is used to combine the standard mineral soil properties 367 with soil organic properties using weighted arithmetic or geometric averages, depending on the 368 parameter (Table 1). An example of this method is shown in Figure 2 for soil porosity, soil 369 saturated hydraulic conductivity and soil heat capacity.  The seasonal and total biases in Table 2 verify this result and show the same behavior for winter 444 soil temperature, although it is difficult to see visually from Figure 5. For these three variables, the 445 simulated time variability is also improved from CTL to SNL to CPT as shown by the other 446 seasonal and total scores (c-rmse and r 2 ) in Table 2. Finally, the new spectral albedo scheme (from 447 CPT to NEW) has a drastic impact on the snowpack simulation in spring. As shown by Figure 5  448 and Table 2, the new spectral albedos clearly improve the simulation of other variables during this 449 period. They induce a sharp springtime snowmelt with a strong decrease in snow depth and SWE. 450 The snow insulation during spring is thus less important and allows the soil surface to warm up 451 faster. As a result, the model is capable of reproducing the strong soil warming observed in April 452 ( Figure 5). Not surprisingly, the soil temperature skill scores for spring and the whole period are 453 drastically improved although there is a slight degradation in winter. 454 Figure 6 shows daily mean time series of the snow density and temperature profiles 455 averaged over the snow season for each experiment. With only 3 snow layers (CTL), the density 456 distribution is more uniform than using the new snow layering scheme with 12 layers (SNL). The 457 significant densification of the bottom layers in SNL is the main process responsible for the snow 458 depth and SWE improvements observed in Figure 5 and Table 2. In addition, the better 459 representation of the vertical profile of density, that induces less dense and thus more insulating 460 surface snow layers from November to February, permits also to better insulate the bottom snow 461 layer from the atmosphere and results in higher bottom snow and top soil temperature. This 462 explains the skill scores improvement found in winter soil temperature in Table 2. The new snow 463 compaction scheme (CPT) tends to increase the density contrast between the top and the bottom 464 snow layers. The snowpack is also denser than with SNL leading to the strong decrease in snow 465 depth observed in Figure 5 and to the better skill scores in snow depth over each period (Table 2). 466 CPT also results in a small warming at the bottom of the snowpack which slightly heats the 467 soil temperature compared to SNL. Finally, the spectral albedo scheme (NEW) has a limited effect 468 on the snow density profile but results in a slightly colder snowpack than in CPT and even SNL 469 (not shown) due to the large daily winter albedos seen in Figure 5. This is the main reason for the 470 lower winter soil temperatures with NEW than CPT and SNL (Table 2). 471 ERA-I elevation is less than 100m to avoid temperature biases for instance that would be directly 499 due to the low resolution of ERA-I. We also only kept the stations where the number of days with 500 a non zero snow depth measurement over the entire period is superior to 100 days, and that have at 501 least 8 days with snow measurement per year. With this filter, the number of available stations 502 decreases to 158, which remains acceptable. Most stations are located in Russia and Western-503

Simulations over North-Eurasia
Siberia with only a few in Eastern-Siberia (Figure 7). bare ground open areas and/or clearings with regular grass cutting following the WMO standards 565 as mentioned previously, they are compared to snow depths and soil temperature profiles 566 simulated by the ISBA bare soil sub-grid patch alone. This patch exhibits conditions which are 567 closest to those at the corresponding field sites, as is generally the case for ISBA in this kind of 568 comparison [Decharme et al. 2013]. The simulation represented here is the NEW-SOC experiment 569 that seems to capture well the snow depth and soil temperature spatial distributions. For snow 570 depth, the latitudinal gradient is well respected. The lower soil temperature along a southwest-571 northeast transect is also well simulated. 572 The seasonal cycles of daily snow depths and monthly soil temperatures (Figure 8 These improvement in snow depth and soil temperature are confirmed by the spatial 591 distributions of their seasonal skill scores (bias and c-rmse). Figure 9 shows the spatial 592 distributions of snow depth seasonal skill scores (bias and c-rmse) during winter and spring. No 593 clear differences among these simulations appear in winter while the bias and c-rmse of many 594 stations are improved in spring by the new snow scheme. The springtime snow depth is simulated 595 in an acceptable manner by NEW, while CTL exhibits a significant overestimation. This fact is 596 confirmed by total scores given in each of the panels. In winter, regardless of the experiments, 597 ISBA underestimates snow depth measurements at many stations, especially in the Northern and 598 Western parts of the domain (Figure 9). 599 The spatial distribution of soil temperature seasonal skill scores simulated at 20 cm and 160 600 cm depth during winter is given in Figure 10. Regardless of the region, the generalized cold bias 601 found over all stations with CTL is drastically reduced with the new snow scheme and the 602 interannual variability (c-rmse) is largely improved. In summer (Figure 11), as was already shown 603 in Figure 8, NEW-SOC is in better agreement with observations compared to NEW regardless of 604 the soil horizon (lower c-rmse) even if a slight cold bias appears at the subsurface as shown by the 605 negative total bias found at 320cm depth. The NEW experiment overestimates the temperature 606 profile measurements at many stations near the surface, but less-so at a 320 cm depth. So, it seems 607 that the subsurface cooling in the NEW-SOC experiment is too intensive. But in fact at 320 cm 608 depth, the simulated soil temperature in the western part of the domain remains quasi unchanged 609 between NEW-SOC and NEW. The best total scores found on Figures 8 and 11  in the soil approximated via a linear interpolation between the last positive temperature node going 620 down from the surface and the first negative temperature node. As expected from the lower 621 summer soil temperatures with NEW-SOC (Figure 9 and 11), the active layer is shallower. 622 However, this comparison with the limits of different permafrost types does not allow to determine 623 which simulation leads to the most accurate active layer thicknesses. The comparison with the 624 CALM data given in Figure 12 seems to show that NEW-SOC simulates a more accurate spatial 625 distribution of the active layer thickness. This result is confirmed by Figure 13  It must be noted that the large improvement in snow albedo in spring is mainly due to the 653 use of snow age in the diagnostics of the optical diameter of snow (Equation 18). Without this 654 parameterization, the surface albedo is strongly overestimated in winter and, to a lesser extent in 655 spring at the Col de Porte field site, with a larger bias and c-rmse for all variables compared to the 656 new version of ISBA (not shown). The optical diameter of snow strongly controls the near-infrared 657 albedo, while impurities mostly affect the albedo in the visible spectrum [Wiscombe and Warren 658 1981]. This increase of snow optical diameter with time is necessary to represent well the decrease 659 in spectrally integrated albedo with age. However, the increase of snow optical diameter is not 660 only a function of snow density as parameterized by Anderson [1976] in Equation (10) The snow thermal conductivity is computed as a function of snow density 731 following Yen [1981]. It also accounts for vapor transfer in the snow using a simple 732 parameterization from Sun et al. [1999]. This process is especially important at low snow densities 733 and at high altitude. So the snow thermal conductivity, λ sn (W.m -1 .K -1 ), in each layer is given by: 734 where ρ snmin = 50kg.m -3 is the minimum density of snow, ρ mob a reference density of 295kg.m -3 , 745 and the dimensionless constant a mob = 1.25. Secondly, a wind-driven compactionindex, Г w , 746 combining the mobility index and the near surface atmospheric wind speed: 747 where κ v = 1.25 is a dimensionless coefficient for gust diagnosis from average wind speed, and the 749 constants a Г = 2.868 and b Г = 0.085 s.m -1 . A positive value of Г w indicates that wind-driven 750 compaction can occur. Compaction rate from the surface is then propagated to the layers beneath, 751 following an exponential decrease, until it meets a snow layer having a negative wind-driven 752 compactionindex. For each layer, this compaction rate is computed as follows: 753 where π τ (s) is a time constant of one day, and the constants a τ = 10 and b τ = 3.25. 755

APPENDIX C 756
Grid-cell snow fraction 757    corresponding skill scores are given in Table 2. Over all panels, the grey shadow corresponds to 1113 the uncertainty in in-situ measurements as discussed in section 4.1. The observed snow depth 1114 exhibits an accuracy of ±1cm, the soil temperature is measured with a precision of ±1K, while 1115 uncertainties in SWE and surface albedo is near ±10%.