Multi-scale validation of a new soil freezing scheme for a land-surface model with physically-based hydrology

Soil freezing is a major feature of boreal regions with substantial impact on climate. The present paper describes the implementation of the thermal and hydrological effects of soil freezing in the land surface model ORCHIDEE, which includes a physical description of continental hydrology. The new soil freezing scheme is evaluated against analytical solutions and in-situ observations at a variety of scales in order to test its numerical robustness, explore its sensitivity to parameterization choices and confront its performance to field measurements at typical application scales. Our soil freezing model exhibits a low sensitivity to the vertical discretization for spatial steps in the range of a few millimetres to a few centimetres. It is however sensitive to the temperature interval around the freezing point where phase change occurs, which should be 1 °C to 2 °C wide. Furthermore, linear and thermodynamical parameterizations of the liquid water content lead to similar results in terms of water redistribution within the soil and thermal evolution under freezing. Our approach does not allow firm discrimination of the performance of one approach over the other. The new soil freezing scheme considerably improves the representation of runoff and river discharge in regions underlain by permafrost or subject to seasonal freezing. A thermodynamical parameterization of the liquid water content appears more appropriate for an integrated description of the hydrological processes at the scale of the vast Siberian basins. The use of a subgrid variability approach and the representation of wetlands could help capture the features of the Arctic hydrological regime with more accuracy. The modeling of the soil thermal regime is generally improved by the representation of soil freezing processes. In particular, the dynamics of the active layer is captured with more accuracy, which is of crucial importance in the prospect of simulations involving the response of frozen carbon stocks to future warming. A realistic simulation of the snow cover and its thermal properties, as well as the representation of an organic horizon with specific thermal and hydrological characteristics, are confirmed to be a pre-requisite for a realistic modeling of the soil thermal dynamics in the Arctic.


Introduction
Frozen soils occupy 55 % to 60 % of the land surface of the Northern Hemisphere in winter (Zhang et al., 2003) with considerable implications for climate (William and Smith, 1989).
Soil freezing impedes water infiltration and drainage, leading to a modified hydrological regime at catchment's scale (Woo et al., 2000). Arctic rivers provide an example of large scale hydrological implications of soil freezing: the Published by Copernicus Publications on behalf of the European Geosciences Union.
seasonal cycle of freshwater input into the Arctic Ocean is highly modulated by terrestrial freeze-thaw cycles ; this freshwater input is of major importance since it partly controls the Arctic Ocean's salinity, sea-ice formation and finally the global thermohaline circulation (McDonald et al., 1999;Peterson et al., 2002;Aagaard and Carmack, 1989;Arnell 2005). In Eurasia, Serreze et al. (2002) found that the runoff to precipitation ratio is proportional to the extent of permafrost underlying river basins. Generally, watersheds underlain by permafrost have a low subsurface water storage capacity (Kane, 1997), implying low winter discharge and a fast hydrological response.
At smaller scales, freeze-thaw cycles induce lateral and vertical water redistribution as a consequence of cryosuction, patterned ground, talik or thermokarst lakes formation. Those features alter the soil structure and thus its water holding capacity, with potential consequences on water fluxes between the soil and the atmosphere, water availability for plants and the functioning of the plant and soil biota ecosystems (Pitman et al., 1999).
Another consequence of soil freezing is the latent heat release and consumption, which delay the seasonal soil temperature signal (Boike et al., 1988). Frozen soils also exhibit specific thermal characteristics due to the different thermal properties of ice and water, and dissimilarities in water distribution within the soil column in frozen and unfrozen states (e.g. Farouki, 1981).
Arctic and boreal regions are in great part underlain by permafrost and/or subject to seasonal freezing. Their soils contain more than 40 % of the global terrestrial carbon (Tarnocai, 2010), undergoing slow or no decomposition due to cold temperatures. The soil microbiological activity is indeed highly sensitive to temperature, especially in subfreezing states (Nobrega et al., 2007). In permafrost regions, microbial activity mostly takes place in summer in the uppermost, thawed layer of the soil, called the active layer. The proper representation of this layer in land surface models is crucial to capture the amount of organic matter decomposition. Depending, among others, on soil moisture conditions, decomposition within the soil will occur through respiration or methanogenesis, leading to the release into the atmosphere of the greenhouse gases CO 2 or methane, respectively. As freezing-thaw cycles and the occurrence of permafrost strongly modulate both the hydric and thermal states of the soil, their representation is crucial in land-surface models including a representation of the carbon cycle. Those models are used in coupled mode to carry out projections of the future climate using different emissions scenarios (IPCC 2007;Friedlingstein et al., 2006). Simulations results reveal the possibility of a strong positive feedback to global warming from the huge high latitude soil carbon reservoir, as increased active layer thickness and permafrost disappearance enhance microbial activity and carbon release to the atmosphere. A representation of the thermal, hydrological and biogeochemical implications of soil freezing is required to improve the realism of those projections (Pitmann et al., 1999;Quinton et al., 2005;Yi et al., 2006).
Soil freezing and permafrost therefore stand out as a key feature for land surface and global climate modeling. Accordingly, efforts have been recently made to introduce a thermal and hydrological parameterization of soil freezing in land-surface models (Luo et al., 2003). Most of these models now include a physically-based representation of hydrology and soil thermal dynamics (e.g. Slater et al., 1998;Smirnova et al., 2000;Essery et al., 2001;Bonan et al., 2002).
The present paper is dedicated to the description and validation of a numerical, one-dimensional soil freezing scheme designed to be part of the physically-based hydrological scheme of the land-surface and carbon model ORCHIDEE (Organizing Carbon and Hydrology Into Dynamical Ecosys-tEms, Krinner et al., 2005). ORCHIDEE commonly provides surface boundary conditions to the atmospheric model LMDZ, but is also used off-line for diverse applications at scales ranging from point location to global. The new soil freezing scheme therefore needs to be validated at a variety of scales. Special attention is given to parameterization and numerical choices and their limits in the context of the current representation of soil freezing in land surface models.
Section 2 fully describes the soil freezing scheme within the model's framework. In Section 3, the scheme is tested against simple, one-dimensional soil freezing experiments and its sensitivity to parameterization choices is discussed. Finally, simulation results at different scales are compared to field data, which helps diagnose improvements induced by the freezing scheme and define further development prospects.

Soil hydrological and thermal processes in the land-surface model ORCHIDEE
ORCHIDEE is the land surface model part of the fully coupled climate model IPSL-CM4 but can be run off-line, driven by a prescribed atmospheric forcing (e.g. reanalyses or outputs from an atmosphere model). It combines a soilvegetation-atmosphere transfer model with a carbon cycle module computing vertically detailed soil carbon dynamics. Although the implications of soil freezing on the carbon cycle are beyond the scope of this paper, the vertically discretized hydrological and carbon modules of ORCHIDEE should provide a useful tool to investigate these interactions. ORCHIDEE computes all the soil-atmosphere-vegetation relevant energy and water exchange processes in 30 min time steps. It is made out of different routines dedicated to energy balance, interaction with the canopy, soil temperatures, soil moisture content, and routing of water towards the oceans. An extended model description can be found in Krinner et al. (2005) for the main land surface processes and in de The Cryosphere, 6, 407-430, 2012 www.the-cryosphere.net/6/407/2012/ Rosnay (1999) and d 'Orgeval et al. (2008) for the vertically discretized hydrology. Hereafter, we only detail the soil hydrological and thermal parameterizations of the model, which are affected by soil freezing. ORCHIDEE allows to choose between a simple hydrological scheme based on 2 reservoirs following the work of Choisnel, (1977) and a vertically discretized hydrological scheme computing vertical water fluxes at each time step (de Rosnay, 1999;d'Orgeval, 2008). Lateral water fluxes are only allowed from one grid-cell to another and do not affect the soil water content, as explained later in the description of the routing scheme. A parameterization of soil freezing exists in the simple hydrology (Koven et al., 2009;Ringeval et al., 2012); however, the improvements induced by a vertically discretized hydrology on the modeling of land-atmosphere water and energy fluxes (de Rosnay, 1999) advocate for the use of a physically-based hydrology and the subsequent implementation of a soil freezing parameterization.
ORCHIDEE vertically discretized hydrology derives from the model of the Centre for Water Resources Research (Dooge et al., 1997) and is documented in de Rosnay (1999) and de Rosnay et al. (2000). It computes the water balance at different depths within the soil profile. Only vertical water movements induced by gravity and suction are accounted for, while water vapour diffusion and water migration driven by osmotic or thermal gradients are ignored. The evolution of soil moisture is thus represented by the 1-dimensional Richards' equation: with : water suction e.g. absolute value of matric potential (m) θ : volumetric water content (m 3 m −3 ) z: depth axis, pointing towards the surface K: hydraulic conductivity (m s −1 ) S: sink term corresponding to water uptake by roots (s −1 ) Equation (1) is discretized on 11 numerical nodes distributed within the soil, with a finer resolution near the surface where key hydrological processes (infiltration, evaporation) take place. These processes affect and in turn are affected by surface soil moisture content. The distance between the two uppermost numerical nodes is 2 mm, and this spatial step increases as a geometric sequence of ratio 2 with increasing depth. The deepest numerical node for the resolution of the hydrological processes is thus at a depth of 2 m. The numerical scheme relies on implicit finite differences and is unconditionally stable. The bottom boundary condition is gravity drainage. At the top of the soil column, the water flux towards the soil is set to infiltration minus evaporation and modulated by the infiltration capacity and water content of the soil. Matric potential and hydraulic conductivity formulations rely on a Van Genuchten (1980) -Mualem (1976) parameterization: with θ s : saturated water content (m 3 m −3 ) θ r : residual water content (m 3 m −3 ) α: Van Genuchten parameter (m −1 ), related to the inverse of the air entry suction m and n: Van Genuchten parameters related to pore-size distribution. m = 1 − 1 n according to the Mualem model. l: Van Genuchten parameter related to tortuosity. l = 0.5 in the Mualem model. K s : saturated hydraulic conductivity (m s −1 ) The parameters α, m, n, and K s are soil-type dependent. The saturated hydraulic conductivity typically varies over several orders of magnitude from coarse to fine-textured soils (Fig. 1a), with considerable impact on the soil water regime. Three different soil types (coarse, medium, and fine) associated with specific hydraulic parameters are accounted for in ORCHIDEE (Table 1). The soil types repartition is the result of the original Food and Agriculture Organization map (1978) and interpolation work by Zobler (1986). In OR-CHIDEE, the original 5 textural classes used by Zobler (fine, medium-fine, medium, medium-coarse, and coarse) are reduced to 3 textural classes (fine, medium, coarse) with the medium class composed of the medium fine, medium, and medium coarse FAO classes. The hydraulic characteristics of the three ORCHIDEE soil textural classes originate from Carsel and Parrish (1988) for the respective USDA (1994) name (Table 1).
Overland flow and drainage water are routed towards the outflow of the major rivers via a routing module thoroughly described in NgoDuc et al. (2007). Basically, the overland flow is transferred to a "fast" reservoir while drainage fuels a "slow" reservoir. Both reservoirs eventually flow into the downstream grid-cell "stream" reservoir, which represents the rivers. The drainage transfer rate from the upstream "slow" reservoir to the downstream "stream" reservoir is slower than the overland flow transfer rate from the upstream "fast" to the downstream "stream" reservoir. The "stream" reservoir water is eventually routed from one gridcell to another till the mouth of the river is reached. All transfer and routing rates depend on the river length from the upstream grid-cell to the down-stream grid-cell and the height loss over that path. The soil temperature is computed according to the Fourier equation using a finite difference implicit scheme with usually 7 numerical nodes unevenly distributed between 0 and 5.5 m. (Hourdin, 1992). The thermal soil is thus thicker than its hydrological counterpart, a necessary feature when considering that the typical damping depth of the temperature annual cycle is about 3 m . The first thermal layer is 4.3 cm thick and the thickness of each layer is multiplied by 2 as the layers get deeper. This resolution was shown to be adapted to the representation of diurnal, annual, and decadal temperature signals (Hourdin, 1992). The upper boundary condition is the flux equilibrium at the soil surface; the lower boundary condition is a zero thermal flux. Latent heat sources and sinks due to the freezing and melting of soil water are by default not included; thermal advection through water movements is neglected. The soil thermal properties depend on the water content, which is interpolated each 30 min time step from the hydrological module at the 7 thermal nodes.

The new soil freezing scheme
The new soil freezing scheme is designed to represent the latent heat exchanges involved in the freezing and melting of soil water, and the changes in thermal and hydrological properties of the ground induced by soil water phase change. Current numerical soil freezing algorithms in land surface models differ in their representation of those effects. The new parameterizations introduced in ORCHIDEE are hereafter detailed and compared with their concurrent counterparts.
Latent heat is a source or a sink term in the Fourier equation. With the assumptions of no thermal advection and no phase change implying the gas phase, the one-dimensional Fourier equation with latent heat term writes: with C p : volumetric soil heat capacity (J K −1 m −3 ) T : soil temperature (K) K th : thermal conductivity (W m −1 K −1 ) ρ ice : ice density (kg m −3 ) L: latent heat of fusion (J kg −1 ) θ ice : volumetric ice content (m 3 m −3 ) In this equation, the mechanical effects of soil freezing (expansion of the total soil volume) are not accounted for.
During freeze-up, latent heat release delays the freezing front progression. Conversely, latent heat consumption counteracts warming as a frozen soil layer reaches the freezing point. As it systematically opposes the temperature change, The Cryosphere, 6, 407-430, 2012 www.the-cryosphere.net/6/407/2012/ latent heat adds up to inertia, which is the basis of its incorporation into an apparent heat capacity in soil freezing models (Fuchs, 1978). This ploy allows to numerically compute a simple diffusion scheme with no source term (Poutou et al., 2004) and is illustrated by the rewriting of Eq. (4) into Eq. (5): with dθ ice dT ≺ 0. The apparent heat capacity can then be analytically derived from the parameterization of the soil volumetric ice content as a function of temperature (Cox et al., 1999;Smirnova et al., 2000). However, numerical complications occur due to the singularity at T = 0 • C. We elude this difficulty following the work of Poutou et al. (2004) with a phase change linearly spread over a 2 • C temperature interval between 0 • C and −2 • C. This temperature interval will hereafter be referred to as the freezing window T . The model sensitivity to the width of the freezing window will be analyzed in Sect. 3. The apparent heat capacity thus simply writes: where θ ice equals the total water content of the layer since all available water is considered to freeze in the freezing window.
To simplify the energy conservation strategy, the total water content used by the thermal scheme does not change from the freezing onset till the end of the thawing. As soon as a layer is entirely thawed, a temperature correction is applied if the amount of latent energy involved in the thawing of this layer and its preceding freezing do not balance each other as a result of numerical approximations. Thermal diffusion is governed by heat capacity and conductivity. The heat capacity of ice is about two times lower than the heat capacity of liquid water; conversely, ice is about four times as conductive as liquid water. Both effects therefore combine to make a thermal signal propagate more rapidly into frozen soil.
The soil heat conductivity K th is calculated as an interpolation between a dry and a saturated value k dry and k sat , on the basis of soil saturation degree.
With k sat : heat conductivity of ice and/or water saturated soil (W m −1 K −1 ) k dry : heat conductivity of dry soil (W m −1 K −1 ) S: total (frozen and unfrozen) soil saturation degree (m 3 m −3 ).
For frozen soils, this parameterization meets up with the Johansen (1975) parameterization, which however recommends the use of 1+log(S) instead of S as a weighing factor for unfrozen soils. We chose S as a unique weighing factor in our parameterization for consistency with the original parameterization of ORCHIDEE (Krinner et al., 2005). The saturated soil heat conductivity k sat depends on soil water and ice content: with k s , k i and k w : heat conductivities for solid soil, ice and water, respectively.
f l : fraction of the liquid soil water, assumed to vary linearly from 1 to 0 between 0 • C and −2 • C, in the "freezing window".
The dry soil heat conductivity is a model parameter and originates from Pielke (2002); its value can be found in Table 2.
The soil heat capacity C p is computed as the sum of the heat capacities of mineral soil and water: With C iw = f l · C wet + (1 − f l ) · C icy : saturated soil heat capacity (J m −3 K −1 ) C dry : dry soil heat capacity (J m −3 K −1 ) C wet : saturated unfrozen soil heat capacity (J m −3 K −1 ) C icy : saturated frozen soil heat capacity (J m −3 K −1 ) The values of these parameters are also listed in Table 2 and originate from Pielke (2002).
Finally, recent studies  have pointed out that an extension of the soil thermal modeling to depths greater than 30 m was needed to prevent unrealistic heat accumulation in the lowest soil thermal layers over decadal to centennial time scales, driven by the zero-flux bottom boundary condition. Simulations over such time-scales are precisely a crucial target for a land-surface model including a representation of permafrost and the carbon cycle, as explained in Sect. 1. In the new soil freezing scheme, the soil thermal column is therefore deepened to 90 m, while maintaining the geometrical increase of the layers' thicknesses: www.the-cryosphere.net/6/407/2012/ The Cryosphere, 6, 407-430, 2012 this vertical extension requires the use of 4 additional layers (Table 3, Extended depth).
The main hydrological impacts of soil freezing are a considerable, though not total, reduction in infiltration and water movements (Burt and William, 1976), concurrent with a low water storage capacity in permafrost regions (Kane, 1997). Those features lead to very specific hydrological regimes in regions underlain with permafrost or subject to long seasonal freezing. Most land-surface schemes assume that water movements within frozen or partially frozen soils occur through unfrozen films and within an unfrozen porosity. These models often prescribe a reduced hydraulic conductivity for frozen soils but still use the Richards' equation to account for water movements. In the SiB and SiB2 model, Sellers et al. (1996b) and Xue et al. (1996) for instance used a linear function to decrease soil hydraulic conductivity at subfreezing temperatures. Lundin (1990) suggested the use of an exponential impedance factor. Other approaches consider that ice becomes part of the soil matrix, which reduces the porosity and the hydraulic conductivity (Kowalczyk et al., 2006). However, this reduction may be too drastic for large scale applications, where water can infiltrate through specific structures like cracks, dead root passages, or where the soil can be locally unfrozen (Koren et al., 1999).
Our new parameterization of frozen soil hydrological processes relies on the two assumptions that (i) only liquid water can move within a frozen or partially frozen soil and (ii) the hydraulic conductivity in a frozen or partially frozen soil depends only on the liquid water content and the soil properties, with no consideration of a reduced porosity due to the presence of ice. The induced reduction in hydraulic conductivity is thus less severe than in most of the above-cited approaches, which could help represent the ability of water to infiltrate frozen soils at a model grid-cell scale through preferential pathways (Koren et al., 1999). This approach furthermore exploits the already available Van Genuchten parameterization of hydraulic conductivity as a function of water content (Eq. (3) and Fig. 1a). Essery and Cox (2001) model the hydrological properties in the land surface model MOSES at subfreezing temperatures similarly.
For the hydrological module, we developed two ways of diagnosing the liquid water content at subfreezing temperatures. The first one, hereafter referred to as "linear" freezing, assumes a linear increase of the frozen water fraction from 0 to 1 in the freezing window, i.e. when temperature goes down from 0 • C to −2 • C. This approach is coherent with the thermal parameterization described above, which assumes a linear phase change for water over the freezing window. The second approach, hereafter referred to as "thermodynamical" freezing, computes the thermodynamically allowed liquid water content at subfreezing temperatures, based on the balance between the low energy status of adsorbed and capillary liquid water, and the free energy drop induced by phase change (Black and Tice, 1989;Dall'Amico, 2010).
We specify here that these diagnostics of the liquid water content are exclusively performed for the hydrological module and do not interfere with the thermal parameterization, which considers a linear phase change over the freezing window regardless of the approach chosen in the hydrological module.
We hereafter describe the thermodynamical approach in more detail. With the assumption of an imposed pressure on ice, Fuchs et al. (1978) derived: With T fr = 273.15 K: freezing point of water? g: standard gravity (m s −2 ) T : soil temperature (K) Equation (10) equally means that soil water under suction will freeze at temperature T ; and if the subfreezing temperature T is observed, the liquid water content has adjusted to the suction .
Liquid water content and soil matric potential are indeed related at subfreezing temperatures, with a relationship similar to what is observed on the course of drying-wetting experiments (Black and Tice, 1989). This suggests that Eq. (2) can be used for frozen or partially frozen soils. A theoretical explanation often advanced (Dall'Amico, 2010) is the replacement of air in the porous media -whose proportion would increase upon drying -by ice when soil freezes. As the stabilizing capillary interactions differ in magnitude between freezing and drying due to a 2.2 times greater surface tension at the air-water than at the ice-water interface, the use of a factor 2.2 in Eq. (2) is sometimes recommended in freezing-thawing applications (Koopmans and Miller, 1966). As capillary interactions are generally involved at lower suctions than adsorptive processes and affect a greater quantity of water, they explain most of the unfrozen water at temperatures just below freezing, when the effects of liquid water are important . The use of the factor 2.2 thus appears relevant, leading to the following equation to describe the thermodynamically allowed liquid water content at subfreezing temperatures as a result of (2) and (10): The real liquid water content is, however, limited by the water available within the soil: With θ l : liquid water content at a subfreezing temperature (m 3 m −3 ) The Cryosphere, 6, 407-430, 2012 www.the-cryosphere.net/6/407/2012/ θ : thermodynamically allowed liquid water content from Eq. (11) (m 3 m −3 ) θ tot : total water content (m 3 m −3 ).
In both the linear and the thermodynamical approaches, the residual water content (see Table 1) does not freeze. Figure 1b and c respectively display the liquid water content diagnosed as a function of temperature by the linear and thermodynamical approaches for the three soil types represented in ORCHIDEE. Fine textured soils retain more liquid water at subfreezing temperatures due to high capillary forces. By contrast, in coarser soils, the decrease in liquid soil water content as a function of temperature is steeper. The simulation was performed with an initial volumetric water content of 0.33 for all soil types at 280 K. Figure 1c also illustrates the limitation of liquid water content by available moisture in coarse soils, since the coarse soil gets depleted in water by gravity drainage before freezing occurs.
The thermodynamical approach is commonly used in land surface models with minor variations (Koren et al., 1999;Cox et al., 1999;Smirnova et al., 2000;Cherkauer and Lettenmeier, 2003). Results yielded by the linear and the thermodynamical approaches will be compared in Sect. 3 and 4.
Other alternative diagnostics of the soil unfrozen water content include a power or modified power function: With b a site specific parameter (e.g. Osterkamp and Romanovsky, 1997); or an ice content determined by total water content and energy loss at T = T fr (Slater et al., 1998;Takata and Kimoto, 2000;Kowalczyk, 2006). We did not try to implement or test them in ORCHIDEE for different reasons: the site-specific calibration requirements disqualify the power function approach for the purpose of land surface modeling at large scales, while the second approach was not easy to conciliate with the original thermal scheme of OR-CHIDEE.

Methods and data
This section aims at evaluating the ability of the new soil freezing scheme to represent the thermal and hydrological processes involved in soil freezing and thawing, and at determining the range of validity of key numerical parameters.
For this purpose, model outputs are compared to idealized data. By "idealized data", we mean data where the unknowns usually restrict the power of model validation (uncertainties in the atmospheric forcing, uncertainties in the soil and vegetation parameters, errors or error compensations due to processes not represented by the model) are minimized. In those conditions, the numerical performance of the algorithm and the suitability of the numerical choices (spatio-temporal discretization, freezing window) can be examined. The mere ability of a scheme to represent a desired process with a desired degree of accuracy is not straightforward as the performance of numerical algorithms are known to be likely sensitive to implementation choices (e.g. de Rosnay et al., 2000;Dall'Amico et al., 2011). A scheme validated in idealized conditions does not necessarily perform well in real climatological conditions. However, establishing the validity and conditions of validity of a numerical scheme is a preliminary step in the validation process.
Regarding soil thermal processes, the analytical solution of the freezing front progression by Stefan (1890)  parameters and boundary conditions used by the model can be set identical to their counterparts in the analytical solution. The accuracy of the scheme can then be investigated in unbiased conditions, and its sensitivity to numerical choices can be explored. To our knowledge, there is no simple analytical solution to the problem of liquid water movements in frozen or partially frozen conditions as described by Eq. (1). We therefore relied on a laboratory experiment to benefit from an explicit setup and well-measured soil parameters which could be used in our model. The Mizoguchi (1990) experiment described in Sect. 3.3 provides such conditions. This experiment is used as a benchmark for the hydrological parameterizations.
This methodology ensures the stepwise validation of the whole soil-freezing scheme: it first focuses on thermal processes alone (see 3.2), then both thermal and hydrological processes are considered, which is necessary due to their intrication. The one-dimensional simulations involved in this section span very limited time ranges (around 50 h) and are of low numerical cost; they are thus particularly suited for sensitivity tests.

Validation of the numerical thermal scheme against the Stefan solution
One-dimensional phase change problems can analytically be solved (Stefan, 1890) with the assumptions of a linear temperature gradient within the soil, a uniform and constant heat conductivity in the frozen zone, and a steady upper boundary condition (Li and Koike, 2003). Those assumptions are rarely met in real conditions. However, they are easily reproduced in the setup of a numerical experiment. In the analytical Stefan solution (Eq. 14), the progression of the freezing front (z) with time (t) is governed by three parameters: the heat conductivity K th of the frozen soil, the surface temperature difference to the freezing point T s -T fr , and the volumetric water content θ of the soil. The latent heat of fusion L and the density of water ρ also appear in the equation.
ORCHIDEE is tested in different configurations against the Stefan solution. A first control run is realized without the freezing scheme (NOFREEZE). Then, a first set of simulations uses the new soil freezing scheme with different vertical discretizations ranging from a regular spacing of 0.005 m between the numerical nodes over the total thermal soil depth (5.5 m), to the default thermal resolution (Table 3; FREEZE, default res). One of those discretizations is set equal to the default vertical discretization of the hydrological scheme (extended to the total thermal soil depth) and will be referred to as FREEZE, improved res. We did not try coarser-thandefault vertical discretizations as the default discretization is dictated by anterior minimum model requirements (Hourdin, 1992). This set of simulations is designed to evaluate the impact of the vertical discretization on the representation of soil freezing. In a second set of simulations, we use the soil freezing scheme with its default vertical discretization and successively set the freezing window T to the following values: 0.1 • C, 0.5 • C, 1 • C and 2 • C. This second set of simulations aims at evaluating the sensitivity of the soil freezing scheme to the freezing window. For all the simulations performed in this section, the time scale involved does not justify the use of an extended depth for the thermal module (Table 3). Furthermore, we also consider a fixed time step of 30 min for the model iterations, as is currently in use in most land surface models and unlikely to change by a factor of more than 2 in the near future.
To suit the conditions for a comparison, ORCHIDEE is set up in conditions close to the Stefan framework: the soil volumetric water content is artificially set at a constant and uniform value, which also ensures constant and uniform thermal conductivity of the frozen soil. Soil temperature is uniformly initialized at 0 • C, and a step-like temperature surface forcing of −6 • C is applied from t = 0 on. The parameters values used for the comparison are: K th = 1.05 W m −1 K −1 ; θ = 0.19 m 3 m −3 ; T s -T fr = 6 K. No hydrological process is involved in the simulations, which allows the testing of the thermal scheme alone. Our numerical experiment only deviates from the Stefan framework by the use of a nonzero heat capacity in ORCHIDEE, while soil heat capacity is not accounted for in the Stefan solution. With our choice of thermal parameters for ORCHIDEE (see above and Table 2), this deviation causes an overestimation of around 8 % (C p ∼2.6.10 6 J m −3 K −1 , C app ∼32.10 6 J m −3 K −1 ) in the energy involved for the phase and temperature change over the freezing window with respect to the Stefan solution. Therefore, a perfect match between our simulation results and the Stefan solution is not ambitioned. We will get back to the consequences of this overestimation at the end of the section. Figure 2 displays the progression of the freezing front within the soil as given by the Stefan solution (STEFAN) and by the NOFREEZE; FREEZE, default res; and FREEZE, improved res simulations described above. In the NOFREEZE case, the progression of the 0 • C isotherm is represented. Both FREEZE simulations considerably improve the representation of the soil thermal dynamics by slowing down the downward progression of the freezing front as compared to NOFREEZE. They also agree well with the Stefan solution at the numerical nodes (RMS < 0.01 m), which are represented by the dashed lines for the FREEZE, default res case. This result is observed for all the finer vertical discretizations we tested (not shown). However, a net overestimation of the freezing front depth, or equivalently a net cold bias, is obvious at depths in-between numerical nodes (step-like features in Fig. 2). This bias is not induced by the soil freezing scheme itself but by the linear interpolation of temperatures between the numerical nodes, as illustrated Fig. 3 can amount to up to 25 % of the analytical solution when the default resolution is used, and is reduced by the use of a finer resolution due to the reduction in the range of the linear interpolation (FREEZE, improved res). Equivalently, the linear interpolation of a summer temperature profile induces a warm bias and an overestimation of the thawing depth also called active layer. A finer-than-default thermal resolution in the uppermost meters of the soil or a specific post-processing of simulation results may therefore be required for specific applications. However, the use of a coarse resolution does not necessarily affect comparisons with active layer thickness observations as the active layer depth is often diagnosed from a linear interpolation of a measured temperature profile (Brown et al., 2003). We conclude that the default vertical discretization, as well as finer resolutions up to 0.005 m, are suitable for the representation of soil freezing, and that finer resolutions may be selected when more precision is required. In Fig. 4, the results of the second set of simulations (sensitivity to the freezing window T ) are analysed in terms of errors in the modeled energy budget. The energy accounted for is the latent heat released upon the complete freezing of the soil water, and it is compared to its theoretical value inferred from the imposed soil moisture content. Figure 4 reveals that freezing windows of 1 • C and 2 • C lead to lower errors than narrower freezing windows at the two soil depths considered. The freezing window can be the source of two types of errors, respectively leading to an underestimation or overestimation of the modeled latent heat. The first one comes from too thin soil layers undergoing temperature changes of higher magnitude than the freezing window and thus overlooking the phase change. This error is responsible for the latent energy deficit in the uppermost 30 cm of the soil with T = 0.1 • C. The second error results from layers whose temperature lies within the freezing window but which undergo a temperature change exceeding the window, thus producing an excess of latent heat in the model. The latent energy overestimation modeled with T = 0.5 • C in the uppermost 30 and 60 cm of the soil is an illustration of this second source of error. Both errors can compensate over time, as illustrated by the case T = 0.1 • C: the uppermost thin soil layers overlook the phase change, which leads to a latent heat deficit in the 30 first cm of the soil. The second error dominates then over the deeper, thicker layers and the error in the latent heat budget is almost corrected when the uppermost 60cm of the soil are considered. Narrow freezing windows and thin layers enhance the freezing-window induced errors. The freezing window should also be coherent with the physics observed: a 0.1 • C-wide freezing window is too small compared to observations (Black and Tice, 1989); a 2 • C-wide window is more realistic when the soil is coarse. Our tests reveal that a freezing window of 2 • C is also compatible with the default vertical resolution of the model. A 2 • C freezing window will therefore be used for the rest of this study.
Nevertheless, a 12 % underestimation of latent heat energy in the uppermost 30 cm of the soil is modeled with the 2 • C freezing window, which is used for the comparison in Fig. 2. As a concomitant 8 % overestimation of heat capacity is inferred by our numerical setup (see above), the good agreement between our numerical solution and the Stefan solution at the numerical nodes in the [0, 30 cm] depth range results from the compensation of both errors. This underlines the limits of both the validation against the Stefan solution and the accuracy of our numerical thermal scheme.

Test of the freezing scheme against the Mizoguchi experiment
Mizoguchi (1990) performed laboratory experiments of soil freezing designed to monitor the evolution of soil moisture as freezing occurs. Four 20 cm deep soil columns of sand of known properties, with initial uniform water content of 0.33, and in thermal equilibrium at 6.7 • C, are placed at t = 0 under a freezing fluid at −6 • C. Only the tops of the columns are sensitive to this boundary condition: the other columns parts are thermally isolated and impermeable. The experiment consists in measuring the soil water distribution after 12, 24 and 48 h of evolution. An unfrozen soil sample serves as a reference. The Mizoguchi experiments also allow the monitoring of the freezing front progression as it corresponds to the zone the most depleted in water. It hence provides a benchmark for the simulation of temperature and water redistribution as a consequence of freezing, in a simplified context where largescale effects or precipitation inputs do not add complexity. The Mizoguchi data were exploited by Hansson et al. (2004) for the evaluation of a numerical heat transport and water flow model. The details of the experimental setup and the hydrological parameters values can be found in this publication.
We created an adapted climatological forcing to test the new soil freezing model against these data. Shortwave radiations were set to zero, incident longwave radiations were chosen as emitted from a blackbody at −6 • C. Wind speed was adapted according to the sensible heat flux transfer coefficient mentioned by Hansson et al., (2004). The model was also configured to suit experimental conditions: the bottom boundary condition was set to zero drainage; the hydrological soil column was limited to 20 cm; the default hydrological vertical discretization was used; the thermal vertical discretization was refined by a factor two, which provided a finer resolution over the 20 cm. Figure 5 compares the results of our soil freezing scheme to the experimental data in terms of freezing front progression (Fig. 5a) and vertical soil moisture distribution (Fig. 5b  and c). Both linear and thermodynamical freezing were tested. A control simulation without the freezing scheme was also performed (not shown), which led to a very slight (0.04/20 cm) vertical soil moisture gradient after 48 h of simulation as a result of hydro-gravitational equilibrium.
The modeled and observed progressions of the freezing front (Fig. 5a) agree well with an error less than 6 % at the 3 time-steps where data are available. This confirms the performance of the thermal scheme. We underline that due to parameterization choices, the freezing front progression modeled with the linear and the thermodynamical freezing do not differ, so that only one model output is plotted on Fig. 5a. The limitations implied by this choice are discussed at the end of this section.
Both the thermodynamical and the linear freezing simulate cryosuction with an amplitude similar to the experimentally observed process (Fig. 5b and c). However, the profiles somehow differ, linear freezing allowing cryosuction to develop deeper within the soil. This can be explained by a less drastic reduction of the liquid water content in the linear freezing when temperature drops below the freezing point ( Fig. 1b and c). On the opposite, cryosuction as modeled by the thermodynamical freezing involves greater liquid soil moisture gradients, which results in water movements of stronger magnitude. These simulations alone do not allow to discriminate the performance of one parameterization over the other.  To our knowledge, validations of the soil freezing hydrology of land surface models against cryosuction data are very scarce (e.g. to some extent Koren et al., 1999). The vertical water redistribution resulting from this process impacts the soil thermal properties and thus the frozen soil thermal dynamics, but the parameterization choices we made do not allow to represent this effect, as the soil moisture used by the thermal scheme remains constant at subfreezing temperature to make energy budget calculations easier. It is also the reason why both soil freezing parameterizations model the same freezing front progression on Fig. 5a. Furthermore, the freezing-induced vertical water redistribution is not expected to have strong implications after soil thawing: in most regions subject to freezing, the uppermost soil experiences saturated conditions in spring as a consequence of snowmelt and/or precipitations infiltration, an effect which is also readily represented by models. This may explain the lack of specific validation attempts of land surface hydrological schemes against cryosuction data. Such a validation however appeared to us meaningful to ascertain the model's physical realism.

Validation at the plot-scale at Valdai
In this section, we use the continuous 18 yrs of atmospheric forcing and hydrological data of the Valdai water balance research station (57.6 • N, 33.1 • E) compiled by Fedorov (1977), Vinnikov et al. (1996) and Schlosser et al. (1997) to evaluate the performance of ORCHIDEE in a region subject to strong seasonal freezing but not underlain by permafrost. These data were extensively used in the PILPS 2d intercomparison project (Schlosser et al., 2000), which provides interesting information about biases in the data and performance of other land surface models. The long-term hydrological measurements are related to the Usadievskiy experimental catchment, whose 0.36 km 2 areal extent is covered with a grassland meadow. The atmospheric forcing data originate from a grassland plot near the catchment; they were initially sampled at 3 h intervals but we used their 30 min interpolation compiled within the frame of PILPS 2d. Incoming longwave radiation was calculated based on the Idso (1981) algorithm. The observed soil parameters for the catchment are extensively described in Schlosser et al. (1997 and. Accordingly, our simulations were performed with a medium soil of rather high hydraulic conductivity (1728 mm d −1 ) and water holding capacity (401 mm m −1 ). The extended vertical discretization was used for the thermal module. Figure 6 compares the mean annual cycles of soil temperature, runoff and soil moisture data over 18 yrs with OR-CHIDEE simulations in three different configurations: without the soil freezing model (NOFREEZE), and with the soil freezing model using the linear parameterization (FREEZE, linear)  The new soil freezing model improves the representation of those three variables. The representation of phase change partially corrects the cold bias of ORCHIDEE in winter and totally offsets its warm bias in summer. The choice between a thermodynamical or a linear parameterization of the liquid moisture content does not impact the modeled soil temperature at 20 cm. As each parameterization leads to a slightly different modeled water content for the 20 first cm of the soil, this result means that the induced soil thermal conductivity and latent energy differences are of minor thermal impact. The remaining winter cold bias possibly originates from (i) the underestimation of the snow cover depth in some winters, as assessed from comparisons to in-situ observations (not shown), (ii) a misevaluation of the thermal parameters of the soil, and/or (iii) the use of a constant, uniform and rather high (0.2 W m K −1 ) snow conductivity. Summer evaporation (and latent heat exchange over the whole year) is marginally impacted by the introduction of soil freezing. The summer soil cooling modeled in the FREEZE-simulations thus originates from a carry-over effect of latent heat consumption during spring thaw in late April. This summer cooling affects the ground below the surface but does not impact the surface temperature itself, which responds to the atmospheric forcing (temperature, radiation).
The timing and amplitude of the runoff spring peak are greatly improved by the soil freezing scheme: in the FREEZE-simulations, the reduced hydraulic conductivity of the frozen soil impedes melt water infiltration, and overland flow is generated when the snow melts in April. When the soil freezing module is not used (NOFREEZE), the spring melt water infiltrates into the soil, leading to a soil recharge visible in the 20 first centimeters of the soil (Fig. 6c). A spring soil recharge of weaker magnitude is recorded in the data, which are averaged from multiple soundings across the catchment (Fig. 6c., DATA). It probably results from macroscale infiltration pathways still active at subfreezing temperatures and from heterogeneous soil freezing conditions at the scale of the catchment, driven by topography and land cover factors. The soil freezing scheme with linear parameterization is able to reproduce a weak spring soil recharge while the thermodynamical parameterization leads to no recharge at all. In the linear parameterization, the reduction of hydraulic conductivity as a function of subfreezing temperature is less drastic and part of the melt water can infiltrate when the temperature is close to the freezing point. Considering infiltration, the linear parameterization therefore appears more suited for an integrated description of catchment-scale processes and heterogeneity. However, the spring soil recharge modeled by this parameterization still occurs one month later than in data, which is in agreement with the delayed spring soil warming illustrated in Fig. 6a. Both FREEZE-parameterizations also capture the autumn runoff peak very well. This peak occurs as a response of a saturated and/or frozen soil to autumn precipitations.
The soil freezing scheme improves the ability of the model to represent the annual amplitude of soil moisture variations. The uppermost soil depletion in water over summer is more marked when the soil freezing module is used, because the late thawing of the deeper, quite wet soil layers in June enhances their hydraulic conductivity and thus the drainage The Cryosphere, 6, 407-430, 2012 www.the-cryosphere.net/6/407/2012/ towards the deeper soil. When freezing is not accounted for, the in-depth soil water profile is closer to hydrodynamic equilibrium and drier at this time of the year; therefore such effects do not occur. In late autumn (mid-November, December), the magnitude of the observed uppermost soil moisture increase is reproduced when the soil freezing model is used, as a cumulated response to autumn precipitation and first freezing events generating cryosuction. The modeled maximum summer depletion in water occurs in August, one month later than in the observed records. This is correlated to a model bias in evapotranspiration, whose summer maximum also displays a one-month lag from the observations (not shown). On average, the soil is wetter than in the observations when the linear parameterization is used. However, this difference is not large (<2.5 % of the average moisture content) and may be due to site-specific characteristics. The reproduction of this phenomenon at other cold-region sites should be verified before we could use it to establish a firm conclusion on the better performance of the thermodynamical parameterization over the linear one.
For soil temperature, runoff and soil moisture, the soil freezing module produces an interannual variability of similar amplitude to the data. The modeled winter soil temperature minima are underestimated due to years where the modeled snow cover is underestimated. The modeled spring runoff minima are also lower than data minima due to years when the timing of the runoff spring peak is not captured accurately by the model. The modeled soil moisture variability is greater in spring and autumn than in data: intricated key hydrological processes occur at these times of the year (freezing, thawing, soil subsurface saturation due to precipitation or snow melt). They exhibit a high spatial heterogeneity and are difficult to capture with a land surface model. Statistics of the modeled and observed interannual variability of soil temperature, moisture and runoff can be found in Table 4. Only the statistics of the thermodynamical parameterization are shown since the linear parameterization leads to comparable results. The modeled interannual variability is improved by the use of the soil freezing module in terms of amplitudes and temporal correlation, except for the soil moisture. However, this last parameter exhibits a rather low interannual variability.
As a conclusion of the Valdai plot-scale model evaluation, the soil freezing scheme noticeably improves the modeled soil thermal and hydrological dynamics at the annual and decadal time-scales. The linear and the thermodynamical parameterizations lead to similar results, with a slightly better representation of spatially integrated infiltration processes in the linear parameterization. We insist that strong conclusions on the soil freezing parameterization can not be drawn from the investigation of one single site, especially regarding hydrology and in link with the high spatial variability of soil hydraulic parameters. Furthermore, the Valdai site may be to some extent representative of regions undergoing sea-sonal freezing. However, it is not located in a permafrost region and therefore limits the extent of our validation to nonpermafrost processes. Further comparisons of model outputs to a full suite of soil temperature and moisture data at cold region and permafrost sites (e.g. described in Westermann et al., 2009;Langer et al., 2011) are required and planned.

Validation across northern Eurasia against soil temperature, active layer and river discharge measurements
We chose to evaluate the new soil freezing model against soil temperature, active layer thickness and river discharge data at the continental scale. Three reasons govern the choice of these variables: (i) they are likely to carry the signature of soil freezing processes; (ii) long-term records exist for them at high latitudes, with an acceptable spatial sampling or spatial representativness (see later in this section); (iii) they are of crucial interest for climate modeling, especially in the prospect of future climate projections. As explained in Sect. 1, the decomposition of organic matter strongly depends on soil temperature. In high latitude soils, it mainly occurs in summer within the active layer, which therefore acts as a key control variable of the high northern latitude carbon balance, with implications for future climate projections (Zimov et al., 2006). Frozen soil processes lead to noticeable changes in the soil moisture regime (see Sect. 4.1) but in-situ soil moisture observations are very scarce, especially in remote high latitude areas. Furthermore the spatial representativness of such measurements is limited (Georgakakos and Baumer, 1996). Conversely, river discharge has been monitored for a long time, especially in the former USSR. This data carries information on the partition between infiltration and runoff, but is spatially integrated at the basinscale. It nevertheless constitutes a crucial climate variable which models should try to represent accurately, since both the amount and timing of freshwater inflow to the ocean systems are important to ocean circulation, salinity and sea ice dynamics (e.g. Peterson et al., 2002). Our study area is northern Eurasia, ranging from 30 • E to 180 • E in longitudes, 45 • N to 80 • N in latitudes. Simulations are performed over this area at a 1 • × 1 • resolution using the meteorological forcing by Sheffield et al. (2006) for the period 1984-2000, with a 10-yr model spinup forced by the 1983 climatology. Soil temperatures are initialized with the mean local ERA (reanalyses) air temperatures over the period . This ensures an initial thermal state as close as possible to present day steady-state conditions and minimizes spinup time requirements for the thermal module. The model is run in three different configurations: without the new soil freezing module (NOFREEZE), with the soil freezing module and the linear parameterization (FREEZE, linear) or the thermodynamical parameterization (FREEZE, thermodynamical). The spatially explicit soil parameters are described in Sect. 1. As the linear and thermodynamical www.the-cryosphere.net/6/407/2012/ The Cryosphere, 6, 407-430, 2012 Table 4. Standard deviation σ and correlation coefficient r between model (FREEZE, thermodynamical) and data for the 20 cm soil temperature, runoff and 20 cm soil moisture over the 18 yrs of simulations and data available at Valdai. For each time period, the statistics are computed using the modeled or observed value averaged over the time period. The statistics improved by the use of the soil freezing module are highlighted.

Soil temperatures
The comparison of soil temperatures simulated by the model with and without the new soil freezing scheme (Fig. 7) highlights the specific signature of the latent heat effects associated with soil freezing. The spring cooling due to latent heat consumption as the soil thaws is visible nearly all over Siberia; the soil thawing occurs later (summer) in the areas with the deepest snow cover (North Western Siberia) as solar radiation energy is first consumed by snow melt. This negative temperature anomaly carries over into the summer. Its magnitude seems less pronounced over the north-eastern coast of Siberia because "summer" encompasses the month of September where the freezing back of the soil has already occurred in these regions. The soil freezing back is responsible for the autumn warming which first affects the coolest, northeastern areas (summer and autumn) and then progresses in a south-western direction along the thermal gradient (winter). In the coldest, north-eastern regions, the winter soil thermal anomaly inverses due to an opposing mechanism: as ice is thermally more conductive than water (see Sect. 2), frozen soils are more conductive than unfrozen ones. In regions with shallow snow cover and mean annual ground temperatures well below the freezing point, at the extent of North-Eastern Siberia, the latter effect dominates over the year and leads to a negative annual temperature anomaly upon the introduction of soil freezing. The same phenomenon is observed in regions with very shallow winter snow cover (Gobi desert for instance), where the poor winter thermal insulation helps the cold wave penetrate faster and deeper within the frozen soil. Over regions with milder winter temperatures or thicker snow cover, the warming of winter soil temperature induced by latent heat effects dominates over the year and leads to increased mean annual ground temperatures. Due to the impact of soil freezing on the soil thermal conductivity (and probably also to hydrological feedbacks), soil freezing hence induces an annual thermal effect although the latent energies involved in freezing and thawing balance out on an annual basis. In-depth soil temperature is monitored at high northern latitudes as part of the Circumpolar Active Layer Monitoring (CALM, Brown et al., 2003) and the Thermal State of the Permafrost (TSP, IPA-SCIDC, 2010). The first dataset consists in measurements of the active layer thickness at CALM stations. The second dataset entails mean annual ground temperatures measured at different depths at the TSP stations. In addition, the Historical Russian Soil Temperature records (HRST, Zhang et al., 2001a) provide a historical perspective on the thermal state of the study area.
We compared the soil temperatures simulated by OR-CHIDEE with and without the soil freezing scheme to HRST records at stations spread over our study area (Fig. 8, top). More precisely, we compared the mean monthly modeled soil temperatures averaged over the decade 1984-1994 to the available data averaged over the same period of time. In the representation of our results, we discriminated the sites where snow is properly modeled from those where it is either underestimated or overestimated. The insulative properties of snow and their crucial impact on the soil thermal regime (e.g. Nicholson and Granberg 1973;Zhang et al., 2005) have indeed long been highlighted by literature, percolating towards the climate (e.g. Luo et al., 2003;Essery et al., 2009) and permafrost (e.g. Kudryavtsev et al., 1974;Zhang et al., 2001b;Stieglitz et al., 2003) modeling communities. Here, the assessment of snow underestimation and the overestimation or correct simulation is based on the comparison of the modeled monthly snow depths averaged over 1984-1994 to the climatological snow depth dataset by Foster and Davy (1988). Generally, the model exhibits a strong cold bias which is reduced but not eliminated by the use of the new soil freezing scheme. We therefore suspect that part of the bias originates from other causes, possibly (i) the misrepresentation of snow (all-year effect) and (ii) the non representation of organic matter (winter effect). These effects will be discussed later in this section.  [1984][1985][1986][1987][1988][1989][1990][1991][1992][1993][1994][1995][1996][1997][1998][1999][2000] period. For each season, the averaged months are mentioned by their initial in brackets. Bottom: annual average over the (1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) period.
The model bias is notably reduced by the new soil freezing scheme in autumn, winter and spring. The latent heat released by the freezing of soil water in autumn warms the soil. This warming can endure over winter and spring, especially if the snow cover, which is a good insulator, is thick. In summer, the soil freezing scheme slightly enhances the cold bias of the model: latent heat consumed by the late thawing of the soils induces a cooling, and the strong model bias in summer proceeds from another mechanism than the representation of soil freezing. The discrimination between the different lev-els of accuracy of the representation of snow in the model (correct, underestimated or overestimated) help point out the role of snow in the thermal regime of the soil in autumn, winter and spring, marginally in summer (Fig. 8): in winter, the model-to-data RMS is reduced by 5.1 K at the sites with correct snow; it is reduced only by 1.7 K at the sites where snow is underestimated. A similar difference in RMS is observed in spring, summer and autumn, though with a reduced amplitude. At those poorly insulated sites, the freezinginduced warming of the soil does not endure over winter as  [1984][1985][1986][1987][1988][1989][1990][1991][1992][1993][1994]. Colors refer to the model with and without the soil new soil freezing scheme; symbols discriminate between the sites where the snow depth is either properly represented by the model within a ±20 cm range (correct snow), or underestimated by the model (snow underestimation), or overestimated by the model (snow overestimation) when compared to the climatological dataset by Foster and Davy (1988).
the uppermost soil is strongly influenced by surface temperatures and its thermal conductivity is enhanced (thermal conductivity of ice is higher than of liquid water). The prominent role of snow in the winter soil thermal regime is confirmed by a forced-snow experiment carried out at Iakutsk, Russia with the new soil freezing model. A 10-yr (1984-1994) simulation is performed using the default snow of the model. In a second simulation, the model is forced by the climatological snow depth extracted from the Foster and Davy (1988)   period). Figure 9 illustrates the changes in the soil thermal regime induced by the use of the climatological snow and compares the modeled soil temperatures to HRST records available during the simulation period. When the climatological snow insulates the soil, the cold bias is strongly reduced at all depths in winter and summer. Although this forced snow experiment is not free of uncertainties (because the forced-snow simulation does not use monthly snow observations over the simulation period but a monthly climatological snow), it still suggests that the poor representation of the Eurasian snow cover could be a major cause of the cold bias affecting the modeled soil temperatures. This was already hinted at by other modeling groups Dankers et al., 2011). The use of a constant and uniform, relatively high thermal conductivity for snow (0.2 W m K −1 ) is another possible contributor to this bias: a high thermal snow conductivity may be adapted for a dense tundra snowpack, but the thermal conductivity of taiga snow is known to be far weaker, with typical values of 0.06 W m K −1 (Sturm and Johnsson 1992). An overestimation of the snow thermal conductivity may explain the increase in the model cold bias from autumn to winter, both with and without the freezing scheme (Fig. 8).
Finally, recent studies report winter soil temperature increases up to 5 • C in boreal regions upon the introduction of an organic horizon into a land surface model (Koven et al., 2009;Rinke et al., 2008;Lawrence and Slater 2008). Organic matter is known for its insulative properties and its specific hydraulic characteristics; it is as well a dominant feature of the Arctic ecosystems (Beringer et al., 2001). Its omission in ORCHIDEE may be a supplementary reason for the model's winter cold bias. However, the abovementioned studies also report significant summer cooling when organic matter is introduced, which would further degrade our model's performance in summer.
Regarding this analysis, we emphasize that the comparison of model gridded outputs to point measurements, as performed here at stations of the HRST network, is of limited relevance due to the extreme spatial variability of soil temperature within a model's grid cell, especially in areas  with complex topography (e.g. Leung and Ghan, 1998). As an illustration, the mean annual ground temperature recorded at 2 TSP stations located ∼15 km apart from each other on the Yamal peninsula differed by 3 K at 1.5 m depth for the year 2008; these two boreholes are located on the same grid cell of the model.

Active layer thickness
Including freezing processes in ORCHIDEE produces an extra-cooling at high northern latitudes, north of 65 • N (Fig. 7). This decreases the modeled active layer thickness over Eurasia, yielding a better agreement with CALM in-situ observations, especially in Eastern Siberia and over the Yamal peninsula (Fig. 10): model-to-data RMS for active layer thickness is reduced from 1.9 m to 0.5 m upon the use of the soil freezing scheme. Our modeled active layer thickness is still overestimated when compared to data, and it is likely that the representation of organic matter in the model would further reduce this model bias, as described in Sect. 4.2.1. and in recent studies (e.g. Nickolsky et al., 2007). Further developments are therefore needed to validate the magnitude of the contribution of soil freezing to the reduction of the model bias in the representation of the active layer. The freezing processes degrade the performance of OR-CHIDEE at most Mongolian sites except in the Altaï range.  The high spatial variability of the active layer in this region, illustrated in Fig. 10, should be kept in mind when evaluating the model's performance; as well as the uncertainties inherent to the comparison of 1 • × 1 • model output to point-scale observation reflecting very local climatological and soil conditions.

River runoff
Soil freezing exhibits a very specific hydrographical signature in regions at least partially underlain by permafrost: in those regions, river discharge is characterized by a very low winter baseflow and a spring peak originating in melt water which does not infiltrate into weakly permeable, frozen soils. The ability of the new soil freezing model to represent these specific dynamics is evaluated by comparisons between modeled and observed hydrographs at the outflow of the three main Siberian rivers Ob, Ienissei and Lena (Fig. 11). The data originates from the Arctic river discharge database R-ArcticNET (Lammers et al., 2001); comparisons are carried out over the 1984-1994 decade, when data are available.
www.the-cryosphere.net/6/407/2012/ The Cryosphere, 6, 407-430, 2012 On the three main Siberian basins, the soil freezing processes similarly impact the modeled hydrographs: the reduction of spring water infiltration within the soil leads to a spring peak of runoff concomitant with the timing of snowmelt. The routing of this overland flow towards the mouth of the rivers, performed by the ORCHIDEE routing module, leads to a spring discharge peak whose timing and magnitude are in agreement with the observed discharge peaks at the outflow of the Lena and the Ienissei. On the contrary, melt water infiltrates within the soil when the physics of soil freezing is not accounted for, and no spring runoff peak is modeled. Drainage occurs at the bottom of the soil, and this subsurface flow sustains the modeled spring discharge peak at the outflow of the rivers (see the dashed lines represented for the Lena basin only in the upper part of Fig. 11). The slower speed of the water flow through subsurface aquifers is responsible for the delay between the spring discharge peaks simulated without and with the soil freezing module, which is also a delay when compared to the observed discharge peaks. We here underline that a spring discharge peak driven by overland flow, as simulated by the soil freezing module, is more reasonable than a drainage-induced spring discharge peak in regions which are partially underlain by permafrost and subject to seasonal freezing. The soil freezing module still does not capture the timing of the spring peak discharge of the Ob river: the vast floodplains of the Ob basin (Ringeval et al., 2010) act as a water reservoir delaying overland flow in its route towards the outflow of the river. These floodplains are not represented in the soil freezing module and may explain the anticipation of the spring discharge peak modeled at the outflow of the Ob. Other possible causes for this wrong timing might be (i) an anticipated timing of the snowmelt over the Ob basin, which is the most temperate region of the study area; and (ii) the water routing scheme, which was not specifically calibrated for Arctic rivers and does not include ice-jam processes. This last reason is however less relevant, as only the Ob discharge timing is miscaptured.
The ability of the soil freezing module to capture the magnitude of the spring discharge peak also varies with the basins and the use of the linear or thermodynamical parameterization. For each river, the linear parameterization leads to a lower spring discharge, because the freezing induced infiltration impedance is less severe in this parameterization and part of the melt water can infiltrate within the soil and be available for either evapotranspiration or subsurface drainage. For the Ob river, both parameterizations of the soil freezing yield a consequent overestimation of the spring discharge magnitude. This bias is probably partly related to the non-representation of floodplains, which in reality foster large evaporation rates and reduce the amount of annual water discharge at the mouth of the river. Floodplains water can also infiltrate the soil later in the year and feed some bottom drainage which sustains the river baseflow, a variable underestimated by our soil freezing model for all basins. For the Lena and Ienissei rivers, the linear parameterization underestimates the magnitude of the spring discharge peak when compared to observations. For the three rivers, the discharge modeled by the linear parameterization also exhibits an unrealistic autumn peak, which originates from bottom drainage at the time of the year where the summer heat wave reaches the bottom of the soil column (2 m) and partly melts the water there, locally increasing the hydraulic conductivity of the soil. This feature is less visible for the Lena basin, which is overall colder than the Ienissei and Ob regions so that soil bottom melting does not occur in autumn; it is also less pronounced with the thermodynamical parameterization for all three basins, because warming closer to the freezing point is required to melt a substantial fraction of soil water in this formulation.
The inability of the soil freezing model to capture the winter baseflow of the main Siberian rivers highlights one of its possible weaknesses. When the soil freezing module is used, autumn rain or melt water hardly infiltrates into the already partially frozen soil and overland flow is produced. In reality, the soil temperatures and thus frozen or unfrozen states exhibit a high spatial variability at the model's gridcell scale (e.g. Leung and Ghan, 1998), and at this scale, part of the water can infiltrate, though with reduced efficiency (e.g. Cherkauer and Lettenmeier, 2003;Niu and Yang, 2006). Taking this subgrid variability into account is likely to sustain a winter, drainage-induced discharge, as simulated when the soil freezing module is not used. The Lena river discharge modeled by the soil freezing module is less affected by this bias: more than 78 % of the Lena basin is underlain by permafrost (Serreze et al., 2002), compared to respectively 36 % and 4 % for the Ienissei and the Ob basins. Accounting for a subgrid variability in the frozen status of the soil is less crucial in the Lena basin since the soil is homogeneously frozen most of the year.
Overall, the new soil freezing model better represents the processes governing the annual dynamics of Siberian rivers. It yields a good agreement between modeled discharge and in-situ data. The thermodynamical parameterization appears more suited for large scale applications. A subgrid variability approach and the representation of floodplains are diagnosed as necessary to capture the annual cycle of the Arctic river discharges with more accuracy.

Conclusion and outlook
A new soil freezing scheme was implemented in the landsurface model ORCHIDEE. It embeds frozen soil processes in a vertically detailed representation of hydrology and the carbon cycle, a crucial feature for carbon-cycle modeling and future climate projections.
The thermal and hydrological behaviors of the new soil freezing scheme are validated at different scales. The scheme thoroughly improves the representation of the soil thermal The Cryosphere, 6, 407-430, 2012 www.the-cryosphere.net/6/407/2012/ dynamics and hydrology at the small and intermediate scales in regions subject to freezing. Over Northern Eurasia, the comparison of present day simulation results to observational data help identify the improvements induced by the representation of soil freezing and the key processes still missed: the soil freezing scheme only partially corrects a winter cold bias in soil temperature, which is partly attributed to an inaccurate representation of snow and the non-representation of organic matter. The representation of the active layer dynamics is more reasonable, yet accounting for organic matter in the model could lead us to reconsider this analysis. The specific features of the Siberian rivers' hydrological regime are captured with an increased accuracy by the new soil freezing scheme. However, the representation of floodplains and the use of a subgrid variability approach appear necessary, especially in regions undergoing seasonal freezing. As a summary, our work highlights the need for development efforts in several directions: -representation of snow (thermal properties and ablation processes) -representation of organic matter (thermal and hydrological effects) -representation of floodplains -representation of a subgrid variability in infiltration These points are currently the focus of new developments in ORCHIDEE.