1The water cycle¶

1.1Hydrological framework: Water fluxes in the soil - plant - atmosphere continuum¶

Ancien text a revoir ou supprimer compte tenu de la section rajoute au tout debut (conceptual framework)

The transpiration sink $s(z,t)$ in equation 5.1 describes the interplay between the transpiration flux, $E_t$, the soil moisture profile, and the root density profile, which is assumed to decrease exponentially with depth, with a decay factor $c_j$ depending on the PFT (Table Y):

Bare soil evaporation $E_g$ is a parallel flux to transpiration, which originates from the entire bare soil column, and from the bare soil fraction of the other soil columns. It is described based on a supply/demand approach:

where the demand is defined by $E_p^*$, the potential evaporation reduced following Milly (1992), and the supply by $Q_{up}$, the maximum amount of water that can be extracted from the soil given the moisture profile (thus at the soil column scale). To prevent from mass conservation violation, $Q_{up}$ is estimated by dummy integrations of equation 5.1 at the end of each time step, in which two boundary conditions are successively tested: firstly a flux condition favoring the demand ($E_g=E_p^*$); then, if the previous case leads to any node with $\theta_i<\theta_r$, a Dirichlet condition ($\theta_1=\theta_r$ at the top node), which strongly limits $E_g$. The resulting water stress factor $\beta m_{G}$ ($E_g=\beta m_{G} E_p$, equation 3.y) will control $E_g$ and the surface energy budget of the next time step. If the total moisture in the top 4 soil layers (assumed to be representative of the litter) is below the wilting point, this stress factor is arbitrarily reduced by a factor 2. Eventually, $E_g$ can proceed at the potential rate of Milly (1992), unless water becomes limiting, i.e. if upward diffusion to the top soil layer cannot provide enough water to sustain the demanded rate.

[AD] This needs to be made consistent with sections 3.1 and 3.4. The case of the bare soil “PFT” needs to be checked, especially for the root decay factor (=5 in the default values). Check if a weighting by the vegetation fractions fvj is required somewhere for transpiration (idem for 1- fvj and bare soil evaporation).

A separate soil water budget is performed in each soil column within a grid-cell. At each time step, the liquid water input to a soil column comprises through-fall and snowmelt, averaged over the contributing PFTs. Given the spatial average of bare soil evaporation and transpiration over the contributing PFTs, the soil hydrology scheme integrates the soil moisture variations over the time step, calculates the output water fluxes (surface runoff and drainage), and diagnoses the water stress factors to calculate the forthcoming evaporative fluxes and feed the STOMATE module (crossref). As detailed in Rosnay (2002) and Campoy et al. (2013), soil water redistribution in the soil implies 1-D vertical water fluxes,

In all cases, $D(\theta)$ and $K(\theta)$ are linearized in the variation range of $\theta$ (between the residual and saturated values, $\theta_r$ and $\theta_s$), by means of 50 piecewise functions, respectively constant and linear in $\theta$ (details in section 1.1.4) This allows a first-order linearization of Equation CITE EQUATION, which can thus be solved by a tridiagonal matrix algorithm. The prognostic variables are the volumetric moisture contents $\theta_i$ of each node $i$ within each soil column. Diagnostic soil moisture layers are also defined, with limits equidistant between two consecutive nodes, so that each layer holds one node. The moisture content of the layers ($W_i$ in mm) is then diagnosed assuming a linear $\theta$ profile between two consecutive nodes.

1.2OK: Canopy interception and througfall¶

As described in Section 1.6.3, the canopy interception flux is calculated for each PFT. It is derived from potential evapotranspiration, moderated by a resistance factor specific to the interception process. This resistance factor incorporates both the aerodynamic and vegetation structural resistances, and is further weighted by the fraction of vegetation that holds intercepted water $f^{inter}$, defined as:

Where $f^{veg}$ is the fraction of vegetation covered by the PFT, $M^{inter}$ is the amount of water intercepted by that PFT per unit of surface (kg m$^{-2}$), and $M^{inter,max}$ is the maximum amount of water that can be intercepted by that PFT per unit of surface (kg m$^{-2}$).

The maximum amount of water intercepted is controlled by the leaf area index $LAI$ (m$^{2}$ m$^{-2}$) and an arbitrary constant $c^{inter} = 0.02$, which transforms $LAI$ into the size of the interception reservoir (unitless), following:

The vegetation canopy is assumed to form a continuous layer, with gaps represented by a supplementary bare soil fraction. For each PFT, the intercepted water storage is first updated by subtracting the interception loss flux over the timestep $F^{ET,inter}$ (Eq. CITE EQUATION Section 1.6.3):

This term represents the loss of water from the interception reservoir due to wet-canopy evaporation. In principle, rainfall first wets the canopy, and only the excess drips to the soil from the intercepted reservoir. In that idealized case, the increase of canopy water storage from a rainfall input $P^{rain}$ would be proportional to the vegetation cover fraction of the PFT: $\delta M^{inter} = f^{veg} P^{rain}$

In reality, however, not all intercepted water remains on the canopy: some evaporates back to the atmosphere (Eq. 150), and some drips to the soil. These processes occur on short time scales (minutes), while the model time step averages precipitation over longer intervals. As a result, the model cannot explicitly resolve the rapid balance between interception, dripping, and evaporation. To account for this limitation, ORCHIDEE prescribes a throughfall fraction $f^{throughfall} = 0.3$ that forces a fraction of rainfall to reach the soil directly. The increase in canopy water storage over the model timestep is then reduced accordingly:

As mentioned previously, canopy storage is limited by a maximum capacity $M^{inter,max}$. When the reservoir exceeds this maximum amount, the canopy water storage is capped and the excess drips into the soil $F^{drip}$:

The input of water to the ground for each PFT $F^{ground}$, therefore includes three contributions: direct throughfall, drip from canopy reservoir overflow, and rainfall onto the unoccupied fraction of the PFT:

Where $f^{veg,max}$ is the maximum fraction of the PFT (unitless).

The volumetric ratio $\frac{M^{inter}}{M^{inter,max}}$ is then a proxy for the actual fraction of the vegetation surface covered with water.

In addition to rainfall, the canopy can also intercept dew (provided it is not freezing). Dew is assumed to behave differently from rainwater: intercepted dew never drips to the soil but always evaporates back to the atmosphere. The fraction of dew that can be intercepted by leaves $g^{dew}$ (unitless), is expressed as a 5-degree polynomial function of $LAI$ CITATION:

This parametrization is only applied for vegetated PFTs with $LAI>1.5$, otherwise $g^{dew}$ is set to 1.

1.3OK: Water status of plants and water transport inside the vegetation¶

The multi-layer soil hydrology scheme presented section 1.7 is partly driven by the vegetation transpiration. This amount of water transpired by the different PFTs (and thus extracted from each soil tiles) is controlled by the rate of aperture of the stomata. This rate of aperture, hereinafter mentioned as stomatal conductance, is regulated by several environmental and physiological factors such as the sensitivity to vapour pressure deficit, leaf sugar concentration and plant hydraulics Damour et al., 2010. This latter stomatal conductance driver, plant hydraulics (and their link to soil moisture), is an essential process to take into account to regulate stomatal aperture and leaves-atmosphere exchanges. It is usually modeled either through an empirical function directly applied to the stomatal conductance formulation Mencuccini et al., 2019Fatichi et al., 2016 or a full representation of the water flow inside the soil-plant-atmosphere continuum which enables the representation of the plant water potentials and a direct control of the stomatal conductance by the leaf water potentials Kennedy et al., 2019Kauwe et al., 2020. In the ORCHIDEE-trunk, both formulations are proposed: i) the historical empirical control of stomatal conductance by soil moisture, described in section 1.3.1; or ii) a plant hydraulics control of stomatal conductance following the implementation of Alléon et al. (submitted), described in section 1.3.3.

1.3.1OK: Empirical dependence¶

In its standard configuration, the stomatal conductance formulation follows Yin & Struik (2009) (more details on the photosynthesis resolution in section 1.2). In this model, stomatal conductance is directly driven by two empirical functions which describe respectively the sensitivity of the stomata to leaf-to-air vapour pressure deficit (approximated here by the air vapour pressure deficit (VPD)) and the sensitivity to soil moisture deficit. For each LAI layer $i$, the stomatal conductance then follows:

where $p^{C_{int}}_i$ and $p^{\Gamma_*}-F^{R_d}_i/k^{g_m}$ represent the inter-cellular $CO_2$ partial pressure and the $CO_2$ compensation point in absence of day respiration ( mol.mol$^{-1}$), $F^A_i$ and $F^{R_d}_i$ represent the assimilation rate and the day respiration ( mol[CO$_2$].m$^{-2}$.s$^{-1}$), $c^{g_0}$ is the residual stomatal conductance (mol[CO$_2$].m$^{-2}$.s$^{-1}$). Finally, G$^{VPD}$ and G$^{hum}$ represent the sensitivities to VPD and soil moisture deficit. They follow respectively:

Where $c^{a_1}$ (unitless, set to 0.85) and $c^{b_1}$ (which equals 0.14 kPa$^{-1}$) drives the sensitivity to $p^{VPD}$, $f^{root}_i$ corresponds to the fraction of root located in soil layer $i$, with $\sum_i f^{root}_i = 1$ $M^{sm}_i$, $M^{smf}$ $M^{smw}$, represent respectively the soil water contents (kg/m$^2$): in layer $i$, at field capacity, and at the wilting point. $M^{sm,nostress}$ is the threshold water content under which stomatal conductance sensitivity $G^{hum}$ starts to decrease linearly and which, with $c^{p_\%}$ (unitless, set to 0.8), puts the start of tranpiration stress at $80\%$ on the segment from $M^{smw}$ to $M^{smf}$. $G^{hum}$ is the sum of the soil moisture stress $g_i^{us}$ across layers, weighted by the root profile.

A non-linear alternative soil water stress can be calculated with the following equation:

ADDED material from Agnes paper

The transpiration sink $S_i$ in Eq. 158 describes the interplay between the transpiration flux, $E_t$, the soil moisture profile, and the root density profile.

The resulting shapes of the transpiration stress factor are illustrated in Figure 6 for the various soil textures considered in ORCHIDEE.

Importantly, all above variables are actually defined at the PFT level, but the index of the PFT was omitted for simplicity. The aggregation at the grid-cell scale is additive, for both the transpiration flux and the total sink term $S_t$, which correspond to the same quantity by construction:

Eventually, in each PFT, the water stress factor $U_s$ is involved in two ways:

• its value from the previous time step defines the root sink term thus the water budget of the current time step;

• once updated at the end of the current time step, based on the corresponding soil moisture, it is used at the beginning of the following time step to calculate the stress function $\beta_3$ controlling the transpiration of the PFT thus the surface energy budget of the following time step.

1.3.2Semi-empirical dependence¶

Describe the Meridja function for humrel

1.3.3OK: Plant hydraulics dependence¶

ORCHIDEE-trunk now embarks a stomatal conductance dependence on leaf water potential. In this configuration, the photosynthesis model of Yin & Struik (2009) is kept but the stomatal conductance formulation is revised and follows Tuzet et al. (2003):

Where $G^{\psi_{leaf}}$ (unitless) represents the sensitivity of stomata to a decrease in leaf water potential ($\psi^{leaf}$ in MPa):

This sigmoïdal function is driven by two PFT parameters: $c^{\psi_{ref}}$ (in MPa), the reference water potential which engenders approximately a 50 % stomatal closure and $c^{s_f}$ (in MPa$^{-1}$), the sensitivity to leaf water potential decrease.

The computation of $\psi^{leaf}$ relies on the representation of water transport from the soil-root interface towards the stomata. In this representation, two water storage pools (for the trunk and the leaves) are defined by capacitances and a set of fixed resistances is defined to model the water flow between each node (i.e root surfaces, trunk-root interface, trunk-trunk storage interface, trunk storage, trunk-leaf storage interface, leaf storage, stomatal cavities) via Ohm’s and Kirchoff’s current laws.

Water storage pools¶

Each water storage pool is defined as in Tuzet et al. (2017) through capacitances ($C^{X,store}$, in m$^3$.MPa$^{-1}$ with X representing either leaf or trunk storage pools) defined according to the storage pool water potential ($\psi^{X,store}$) by the following equation:

With,

Where, $M^{X,max}$ (m$^3$) is the maximum volume that can be stored in the pool (arbitrarily defined as 80 % of the biomass of carbon of the leaves and the sapwood as described in section 1.4.4), $M^{X,res}$ (m$^3$) is the residual volume of water of the pool (arbitrarily defined as 40 % of $M^{X,max}$ after personal communication with A. Tuzet), $c^{\lambda}$ (MPa$^{-1}$) is, according to Tuzet et al. (2017) a parameter depending on the morphological adaptation in different species and on environmental conditions experienced during growth and $c^{\psi_0}$ (MPa) is a reference water potential.

The dynamics of each storage pool rely on a differential equation (Eq. 164) which links the storage compartment water potential ($\psi^{X,store}$) to the flux between the trunk-storage compartment interface and node above ($F^{up}$, in m$^3$.s$^{-1}$), to the resistances above and below the storage compartment ($R^{above}$ and $R^{below}$, in MPa.s.m$^{-3}$) and to the storage compartment resistance ($R^{X,store}$, in MPa.s.m$^{-3}$).

Both storage pool associated differential equations are solved through a predictor-corrector scheme (a two-steps resolution where the state variable at time step $t+1$ is firstly predicted with a classical explicit resolution and then corrected by a second calculation involving the prediction). This resolution permits to avoid instabilities when calculating the new value of the storage water potential. In this resolution, the predictor term $\{\psi^{X,store}\}^{t,predict}$ is calculated and used as follows:

Where, $f(t,\psi^{X,store}_t)$ corresponds to the right-hand side of Eq. 164.

This predictor-corrector scheme enables to approximate the value of the storage pool water potential at the end of the time step which, along with the value at the beginning of the time step, leads to the calculation of the water fluxes coming from or inside the storage pool.

Water absorption by roots¶

At the soil-root interface, the water absorption by roots is computed in two soil horizons (averaged from the multi-layer soil hydrology scheme presented in section 1.8). It can be represented by two different approaches.

In the simpler approach, a dynamic resistance is applied at the soil root interface. Its formulation (Eq. 167 for a soil layer $i$) follows the representation implemented by Bonan et al. (2014) (Section A4.3, Equations A23 to A25).

Where, $L^{roots}_i$ is the total root length per unit volume of soil in the soil horizon (m.m$^{-3}$ - calculated according to root biomass), $\Delta z_i$ is the soil horizon height (in m), $r^{roots}$ is the mean fine roots radius (which equals 0.5 mm) and $k^s_i$ is the hydraulic conductivity in layer $i$ (in mm.day$^{-1}$).

In this framework, the soil horizon water potential and hydraulic conductivity are calculated following either Genuchten (1980) formulation or Campbell (1986) in order to avoid the numerical instabilities which can appear when the soil water content is close to the residual soil water content.

In the more complex framework, the water absorption by roots is represented following Tuzet et al. (2003). In this mechanistic representation, the volume of soil in the horizon is considered cylindrically distributed around a fictive root of length $L_{roots_i}$. In this cylinder of soil, Richard’s equation (Eq. 168) is solved radially which enables to compute the gradient of moisture content near the root.

where $D(\theta)$ (mm$^2$.day$^{-1}$) is the soil diffusivity and $r$, the distance from the root axis. Richard’s equation is solved at a half hourly time step for $n^{muff}$ nodes along the radial axis. At the soil-root interface ($i=0$), the absorbed water flux ($F^{root,up}$ or $F^{root,low}$, in m$^3$.s$^{-1}$, respectively for the upper and lower soil horizons) is applied as a boundary condition. The radial resolution of Richards’ equation is performed using the same resolution framework as the vertical water diffusion presented before (see Section 1.7). This resolution enables the calculation of the water potentials at the root surfaces ($\psi^{rootsurface,up}$ and $\psi^{rootsurface,low}$), by considering a continuity between the soil water potential at the node $i=0$ and the root surface.

Numerical resolution¶

During unstressed periods, the resolution of the hydraulic architecture is purely explicit. At the beginning of the time step, photosynthesis (section 1.2) is solved and $g^s_i$ is the first calculated variable, using Eq. 160 and the leaf water potential from the previous time step. Then, the resolution of the energy budget (sections 1.1 & 1.8) enables the calculation of the transpiration flux, $F^{transp}$ using $g^s_i$. The different water fluxes within the hydraulic system are computed from the top to the bottom imposing $F^{transp}$ flux as a boundary condition at the leaf level. The flux between each storage pool and the nearest node ($F^{trunk,store}$ or $F^{leaf,store}$) is assessed following the predictor-corrector scheme described in the previous section. Using Kirchhoff’s current law and the storage water fluxes, all water fluxes can be determined down to $F^{root-trunk}$, the flux between the trunk and the trunk-root interface, above the two root horizons. At this interface, the two fluxes of water uptake by roots, $F^{root,up}$ and $F^{root,low}$, are assessed from $F^{root-trunk}$ and the root water potentials of the previous time step with Eq. 169.

Based on the two root water fluxes, water absorption by roots is solved. The water contents in the two soil layers are updated as well as the resulting root water potentials $\psi^{rootsurf,up}$ and $\psi^{rootsurf,low}$. The other water potentials of the hydraulic system are then updated, from the bottom to the top, using the previously calculated fluxes and the water potential of the compartment down the flux. Thus, the updated water potential at a given level ($\psi^{upper,level}$) is expressed as:

Where $\psi^{lower,level}$ is the water potential at the level down the flux; and $F^{level}$ and $R^{level}$ are the flux and the resistance between the two levels, respectively.

During drought periods, the non linearity of the stomatal conductance formulation and the strong coupling between stomatal conductance, transpiration and leaf water potential can introduce strong feedbacks between the variables, further leading to numerical instabilities which would justify the use of an iterative resolution scheme (Alléon et al. (submitted), Yang et al. (2019)Sabot et al. (2022)). In ORCHIDEE-trunk, this drawback has been partly solved by implementing an adapted resolution scheme.

A first resolution of the hydraulic architecture model is performed using the explicit framework presented previously. This first calculation is purely diagnostic and enables the calculation of the variables at the current time step. A second resolution of the hydraulic architecture model is then launched. The objective of this second calculation is to approximate the value of the leaf water potential at the next time step in order to avoid a resulting stomatal conductance too far from reality at the next time step. To do this, an approximation of the transpiration flux at the next time step is performed. The transpiration equation can be summarized as:

where $\rho$ is the air density (kg.m$^{-3}$), $\Delta t$ is the time step length, $\Delta q$ is the difference of humidity between the surface and the atmosphere (kg.kg$^{-1}$), $R^{aero} = \frac{1}{\vert \overrightarrow{V}\vert C_d}$ is the aerodynamic resistance (m.s$^{-1}$) and $R^g = 1/\sum_i g^s_i\cdot f^{veg}_i$ is the stomatal resistance (m.s$^{-1}$). This equation can be simplified to:

The transpiration at time step $t$ and at time step $t+1$ can be expressed as:

The objective here is to estimate $F^{transp}_{t+1}$ from $F^{transp}_t$. Looking at the components of the transpiration:

• $A_t$ and $A_{t+1}$ represent the evolution of $\Delta q$ over the day;

• $B_t$ and $B_{t+1}$ represent the evolution of the aerodynamic resistance;

• $C_t$ and $C_{t+1}$ represent the evolution of the components of $g^s_i$ that are not directly linked to $g^{\psi}$.

The following assumptions are made:

• $A_{t+1}$ = $A_t \cdot \delta$;

• $C_{t+1}$ = $C_t \cdot \delta$;

• $B_{t+1}$ = $B_t$;

• $\delta$ can be estimated by a sinusoidal function centered on solar noon.

The assumption is, hence, that in the absence of water stress, transpiration and stomatal conductance will follow a sinusoidal function centered on the solar noon, starting at sunrise and ending at sunset. Following this assumption, $\delta$ can be expressed as the ratio between the values of the sinusoidal function at time step $t$ and $t+1$. The amplitude of the sine function cancels out in the ratio.

Consequently, the ratio between $F^{transp}_{t+1}$ and $F^{transp}_t$ can be estimated as:

This ratio, multiplied by $F^{transp}_t$, permits to have a prediction of $F^{transp}_{t+1}$ only with variables from the previous time step. This approximation enables the second solving of the hydraulic architecture and the estimation of the leaf water potential at the next time step. As the transpiration estimation is relying on the leaf water potential value at the same time step, several iterations are performed in order to make the resolution converge towards a unique leaf water potential value which will be sent to the photosynthesis module at the next time step. One can note that the scheme can be improved by improving the transpiration flux estimation.

1.4OK: Snow and land ice¶

The multi-layer snow scheme described and evaluated in Wang et al. (2013) is based on the ISBA-ES model (Boone & Etchevers (2001)) and has been implemented on the vegetated/biological fraction of the ORCHIDEE grid cell. Its extension to glaciers has been validated over Greenland ice-sheet (Charbit et al. (2024)) and replace the old single-layer scheme described in Chalita and Letreut, 1994. This multi-layer scheme is used to describe the hydro-thermal processes within the snowpack and solve the mass and energy balance equations.

To describe the snowpack evolution, a number of three layers was chosen, for which three prognostic variables (i.e., snow heat content, density and thickness) are calculated. Temperature and liquid water content are diagnosed at each time step and for each layer. The total snow mass and its mean density is then used to derive a snow cover fraction that enters the calculation of the surface albedo and roughness of each grid cell that contains snow, to account for the high albedo and smoothing impacts of snow Wang et al. (2013).

The multi-layer snow scheme describes the energy and mass transfers between the snowpack and the atmosphere and the main internal processes at play such as snow settling, water percolation and refreezing. The calculations are done in the following way: at each time-step, the surface flux entering the top snow layer is computed from the surface energy balance equation. Snowfall updates the thickness, density, heat content, temperature and water content of the snow layers. Compaction is then modeled, modifying layer thicknesses, densities and heat content (equations (11)-(12) in Wang et al. (2013)). Temperature is then updated, from top to bottom, taking into account the new surface temperature and the integration coefficients for the snow numerical scheme, similarly to the soil temperature scheme (ref Hourdin yyyy ?) to stick with the implicit scheme. Snowmelt is then calculated, followed by liquid water transfers and possible refreezing, energy transfers during these phase changes are accounted to update temperature, layer thickness, density and liquid water content. Sublimation is then calculated followed by snow aging, snow cover fraction and finally the surface albedo as presented in section ??

1.4.1OK: Snow energy balance¶

The surface flux entering the top snow layer is computed from the surface energy balance equation in the same way as for soil. At the surface of the snow, the energy budget equation is:

Where $G1$ is the soil heat flux due to heat conduction process; $H1$ is the energy released by rainfall (see Eq. (14) in Boone & Etchevers (2001)). $Frad$, $F1h$, and $LF1$ are the net radiation, sensible heat and latent heat flux respectively ($W\,m^{-2}$); $C_{surf}$ is the surface heat capacity of soil per unit area ($J\,m^{-2}\,K^{-1}$) and is computed as the sum of heat capacities for snow-covered and snow-free surfaces weighted by their respective grid cell fraction.

The heat conduction flux into the soil $G1$ is used to compute the surface temperature $T_{surf}$ of the grid cell at the next time step and provides the limit condition of the surface temperature at the snow-atmosphere interface for the calculation of the snow temperature profile.

Above snow-covered surfaces, when $T_{surf}$ is above the freezing temperature $T_0$ (273.15 K), the energy excess is first used to bring the snow temperature to $T_0$. A surface energy flux $F_{freezing}$ associated with the freezing temperature $T_0$ can be computed using a similar formulation to Eq. 176. The difference between $T_{surf}$ and $T_{freezing}$ is converted in an additional temperature expressed as:

If $T_{add\_snow}$ is greater than $T_0$, the energy excess is used for melting snow, and $G1$ is further set to $G_{freezing}$ for energy conservation. If the new $G1$ value is greater than the total heat content of the snowpack, snow is entirely melted and the excess energy is transferred to the underlying soil. The energy released by snowfall is accounted for in the snowpack scheme to update the snow heat content of the snowpack after a snowfall event.

The sensible $F1h$ and latent heat $LF1$ fluxes computed for each grid cell are given respectively by:

where $\rho_{air}$ is the air density, $T_{surf}$ and $T_{atm}$ are the surface and the 2 m atmospheric temperatures, $Q_{air}$ and $Q_{sat}$ are the air specific humidity at 2 m and the saturated specific humidity at the surface, $L_s$ is the latent heat of sublimation (2.8345 106 J kg-1), U is the wind speed at 10 m and $q_{cdrag}$ is the drag coefficient computed as a function of the ice roughness length ($z0_{ice} = 0.001 m$), following the Monin-Obukhov turbulence similarity theory (Monin & Obukhov (1954)) and the parameterizations of the eddy fluxes proposed by Louis (1979).

1.4.2OK: Snow layering¶

In the initial version of ORCHIDEE explicit snow, the snow was distributed over 3 layers, with the first layer limited to a depth of 5cm and the second to 50cm. This vertical discretisation described in Wang et al. (2013) is still available but we now use a new discretisation with 12 layers. With this new version, the maximum thickness of the first 3 layers is respectively 1, 5 and 15cm, which limits the bias on the surface temperature. An application to the Greenland ice sheet Charbit et al. (2024) has shown that this new 12 layers scheme provides a better simulation of snow evolution than the previous 3 layers model. The definition of snow layer thickness followed the layering scheme proposed by Decharme et al. (2016) for ISBA-ES: In theses equations, $i$ identifies the snow layer and D the snow layer thickness

where the constants $\delta_i$ are the maximum snow depth for a layer $i$ with $\delta_1=0.01\,m$, $\delta_2=0.05\,m$, $\delta_3=0.15\,m$, $\delta_4=0.5\,m$, $\delta_5=1\,m$, $\delta_9=1\,m$, $\delta_{10}=0.5\,m$, $\delta_{11}=0.1\,m$, $\delta_{12}=0.02\,m$. For very thin snowpack (less than 0.1 m), each layer has the same thickness and when the snow depth is above 0.2 m, the first and last layers reach their constant values of 0.01 and 0.02 m. The layer thickness are updated if one of the first two layers or the bottom layer become too thin or too thick:

For example, if the thickness of the top layer becomes lower than 0.005 m or greater than 0.015 m, all the layer thicknesses are recalculated with Eq. 180 and the snow mass and heat content are redistributed accordingly.

1.4.3Ok: Snow compaction¶

When a snowfall occurs, the snow depth increases, but the depth decreases with compaction. This is why compaction is a particularly important process for the evolution of the density and depth of snow in each layer. Two processes are responsible for compaction: the weight of the overlying snow layers and the metamorphism of the snow. These two processes are taken into account in the following equation derived from Anderson (1976):

The first term of the right-hand side of equation 182 represents the compaction due to snow load with $\sigma_i$ (Pa) being the pressure of the overlying snow and $\eta_i$ (Pa.s) the snow viscosity. The viscosity is an exponential function of snow temperature and density (equation 183, Mallor (1964) and Kojima (1967)):

with $\eta_0$=3.7 10$^{7}$ Pa s, $a_{\eta}$=8.1 10$^{-2}$ K$^{-1}$ and $b_{\eta}$=1.8 10$^{-2}$ m$^3$ kg$^{-1}$.

The second term $\xi$ of the right-hand side of equation 182 represents the effect of metamorphism:

with a$_{\xi}$ = 2.8 10$^{-6}$ s$^{-1}$, b$_{\xi}$ = 4.2 10$^{-2}$ K$^{-1}$, c$_{\xi}$ = 460 kg m$^{-3}$ and $\rho_{\xi}$ = 150 kg m$^{-3}$. In the model, the density cannot exceed a threshold fixed at 750 kg m$^{-3}$. Compaction does not affect the total mass and the heat content of the snowpack but changes the layer thicknesses, this is the reason why the snow heat within the layers must therefore be updated.

1.4.3.1OK: Heat content¶

For processes linked to phase change, we use the equation of the heat content of snow (equation 185). This allows to determine whether the snow is cold (dry) or warm (wet). The heat content of each snow layer is computed using the following equation:

where L$_f$ is the latent heat of fusion and T$_f$ the triple-point temperature for water. H$_i$ is the total heat content, c$_i$ the snow heat capacity (J.m-2 K-1) and W$_{liq,i}$ liquid content for the i$^{th}$ layer. According to equation 185, heat content is used to diagnose the snow temperature assuming that there is no liquid water in the snow layer. If the calculated snow temperature exceeds the freezing point, snow temperature is set to T$_f$ and the liquid water content is diagnosed from equation 185.

1.4.3.2OK: Melt and refreezing¶

When melting occurs at the surface, the liquid water remains in the first layer if the maximum liquid water threshold has not been exceeded. When the amount of liquid water exceeds this threshold, the water percolates into the lower layer. The water can then either refreeze or be added to the existing liquid water, depending on the temperature of the snow layer. The maximum liquid water holding capacity is taken as a function of the snow layer density following Anderson (1976).

with $r_{min}$=0.03, $r_{max}$=0.10 and $\rho_t$=200 kg m$^{-3}$.

When the liquid water reaches the bottom layer and exceeds the total liquid water content, it is then considered as runoff.

1.4.3.3OK: Heat conduction¶

The heat conduction from the surface to the bottom of the snowpack is described by a vertical diffusion equation relating the temporal evolution of the snow temperature in the snowpack at a depth z and the divergence of the snow heat flux F$_c$ and is solved using an implicit numerical scheme.

with $\Lambda$ the snow thermal conductivity (W m$^{-1}$ K$^{-1}$). At the snow-atmosphere interface, the boundary condition is given by the energy balance equation 176.

Along with the thermal gradient, a water vapor diffusive flux takes place from the warmer to the colder parts of the snowpack and sublimation or condensation may occur in the pore spaces depending on the water vapor saturation pressure. This process is particularly significant in the Arctic because of strong temperature gradients between soils and atmosphere and is in great part responsible for snow metamorphism. While it is explicitly accounted for in detailed snow models, in Explicit Snow, the effect of water vapor diffusion and phase changes is parameterized through the thermal conductivity (Sun et al. (1999)). The thermal conductivity $\Lambda$ is expressed as the sum of empirical formulations for snow thermal conductivity $\Lambda_{cond}$ (eq. 189) and thermal conductivity from vapor transport $\Lambda_{vap}$ (eq. 190) with:

where a$_{\lambda c}$ = 0.02 W m$^{-1}$ K$^{-1}$, b$_{\lambda c}$ = 2.5 10$^{-6}$ W m$^{-1}$ K$^{-1}$ (Anderson (1976)), a$_{\lambda v}$ = -0.06023 W m$^{-1}$ K$^{-1}$, b$_{\lambda v}$ = -2.5425 W m$^{-1}$ and c$_{\lambda v}$ = -289.99 K (Yen (1981)). P is the atmospheric pressure in hPa and P$_0$ = 1000 hPa.

1.4.3.4OK: Ice layers¶

In previous versions of ORCHIDEE, ice-covered surfaces were simulated using the single layer snow scheme described in Chalita & Le Treut (1994). Now, the multi layer snow model is used for all surface types. Furthermore, in order to be able to simulate the surface mass balance on glaciers and ice sheets, the treatment of ice-covered areas has been modified (Charbit et al. (2024)). To take account of the presence of ice, we have added an ice module between the ground and the snow. This ice reservoir is made up of 8 layers of fixed thickness (0.01, 0.05, 0.15, 0.5, 1, 5, 10 and 50 m). A finer vertical spacing is imposed for the upper layers to better resolve heat conduction at the snow-ice or atmosphere-ice interface. The large thickness of the bottom layer allows it to have an almost constant temperature throughout the year as it has been observed at a few tens of meters depth (Patterson (1994)). The temperature in ice is obtained using the same numerical scheme as for the snow. The formulation of heat capacity $C_{ice}$ and thermal conductivity $\lambda_{ice}$ of the ice are coming from Yen (1981).

where $T_{ice}$ is the ice temperature, $a_{c_ice}$=2115.3 J K$^{-1}$ kg$^{-1}$, $b_{c_ice}$ = 7.79293 J K$^{-2}$ kg$^{-1}$, $a_{\lambda_ice}$= 6.727 W m$^{-1}$ K$^{-1}$ and $b_{\lambda_ice}$= -0.041 K$^{-1}$. Ice is considered as an impermeable medium, hence liquid water coming from melting ice is considered to runoff instantaneously with no possibility of refreezing.

1.4.3.5OK: Snow cover fraction¶

Snow cover fraction (SCF) is parameterized according to snow average thickness and density following Niu and Yang (2007) with revised values for minimal density of fresh snow $\rho_{min}$ (prescribed now to 50 kg/m2), ground roughness length $z_{0g}$ equal to 0.01 m, and a value of 1 for the adjustable parameter “m”.

The thermal properties such as the heat capacity and thermal conductivity are also computed for each snow layer to describe the heat conduction within the snowpack and the temperature profile. (equation (5) in Wang et al. (2013)). Finally, the integration coefficients for the snow thermal scheme are computed, just after the ones for the soil in the same bottom-up direction. This ensures that the multi-layer snow module completely respects the implicit scheme.

SCF is then used to weight the albedo and surface roughness of the grid cell containing glaciers. Snowmelt and runoff are calculated at each time step when surface temperature is above freezing. Soil thermics is calculated by replacing the uppermost soil layers by snow according to snow depth. Heat capacity and thermal conductivity of snow are used for the layers filled by snow.

1.5Soil evaporation¶

$Q_0$ (see XXX) is defined by the difference between infiltration into the soil (section 1.7.10) and soil evaporation, $E_g$.

The latter is a parallel flux to transpiration, which originates from the entire bare soil column, and from the bare soil fraction of the other soil columns. It is calculated using a supply/demand approach, assuming it can proceed at the potential rate (Eq. %s), unless water becomes limiting, i.e. if the upward diffusion to the top soil layer cannot provide enough water to sustain the required potential rate. If this is not the case, the diffusion needs to be recalculated assuming a lower soil evaporation.

*** Explain a bit the correction of Milly *** The potential rate is now calculated using the correction of milly1992potential. The corrected potential evaporation is further called $E^*_{\rm{pot}}$.

The limitation of bare soil evaporation by upward diffusion is estimated based on a dummy integration of Richards redistribution. *** better articlulate with the following

In each soil column $c$, we define the stress factor $\beta_g^c$ aiming at controlling the calculation of bare soil evaporation $E_g^c$, by

The demand is defined by $E^*_{\rm{pot}}$, the potential evaporation reduced following milly1992potential, and the supply by $Q^c_{\rm{up}}$, the maximum amount of water that can be extracted from the soil column given the moisture profile (thus at the soil column scale).

To prevent from mass conservation violation, $Q^c_{\rm{up}}$ is estimated by dummy integrations of the Richards redistribution equation at the end of each time step, in which two boundary conditions are successively tested: