1The energy budget¶

1.1OK: The surface energy budget: principle and main equation¶

ORCHIDEE’s energy budget follows a "big-leaf" approach, in which the surface (of temperature $T^{surf}$ (K) and saturated humidity $q^{surf} = g^{q_{sat}}(T^{surf})$ ($kg.kg^{-1}$)) is considered as a layer of infinitesimal thickness that exchanges with the atmosphere. The energy budget of the surface is based on two main equations. The first corresponds to the radiative budget at the surface:

Where, $F^{R_n}$ (W/m$^2$) is the net radiation at the surface, $F^{LW\downarrow}$ and $F^{SW\downarrow}$ (W/m$^2$) are the downwelling long- and short-wave radiations received respectively from the atmosphere and the sun, and $F^{LW\uparrow}$ and $F^{SW\uparrow}$ (W/m$^2$) are the upwelling long- and short-wave radiations emitted and reflected by the surface.

The second equation represents how the energy available at the surface ($F^{R_n}$) is used. By heating the surface, the absorbed energy is either re-emitted to the atmosphere under the form of convective fluxes, latent and sensible heat fluxes ($F^{LE}$ and $F^{SH}$ respectively - $W/m^2$) or transmitted to the ground through the ground heat flux ($F^G$ - $W/m^2$). The resulting equation enables us to calculate the evolution of the temperature of the surface $T^{surf}$ ($c^{cp}$ being the surface heat capacity $\mathrm{m}^2.\mathrm{s}^{-2}.\mathrm{K}^{-1}$).

At equilibrium, the surface energy budget equation corresponds to:

1.2OK: Radiative transfers¶

As presented in section 1.1, the radiative budget equation allows us to calculate the amount of energy absorbed by and available at the surface. It is composed of net short- and long-wave radiation budgets.

1.2.1OK: Short-wave radiations budget¶

The short-wave radiation budget calculates the amount of solar radiation absorbed by the surface. Incoming short-wave radiation ($F^{SW\downarrow}$) is an input variable of ORCHIDEE. The upwelling short-wave radiation ($F^{SW\uparrow}$) represents the amount of solar radiation that is reflected to the atmosphere, as a function of the albedo ($\alpha$) of the entire grid-cell. Section 1.3 details the computation of the albedo for all sub-grid tiles.

1.2.2OK: Long-wave radiations budget¶

The incoming long-wave radiation emitted by the atmosphere and received by the surface ($F^{LW\downarrow}$) is, as for the incoming solar radiations, an input of ORCHIDEE.

The upwelling long-wave radiation ($F^{LW\uparrow}$) represents the emissions of the Earth surface, which is assumed to behave as a grey-body of emissivity $c^\epsilon = 1 ??$:

Where $c^\sigma = 5.67.10^{-8}$ W.m$^{-2}$.K$^{-4}$ is the Stefan-Boltzmann constant.

1.3Surface albedo calculation¶

1.3.1Grid cell albedo¶

The albedo scheme divides the solar spectrum into two wavelength bands: visible (VIS) and near-infrared (NIR), with a separation at XX nm. For each domain, the overall grid cell albedo result from the combination of the albedo of specific grid cell components. Overall, the surface is divided into three compartments: bare soil with fractional cover $f^{bs}$, vegetated surfaces with fractional cover $f^{veg}$, non-biological surfaces with fraction cover $f^{nobio}$.

For the vegetated cover, the maximum surface occupied by each PFT is divided into two parts for the radiative transfer scheme (see section below): a fraction covered by leaves and a fraction corresponding to bare soil which represents gaps in the canopy or the absence of leaves. These fractions thus evolve over time. For each grid cell the different bare soil fractions are grouped together with the same albedo, i.e. the so-called background albedo. Overall the surface is thus divided into three compartments: bare soil (bs), vegetation cover (veg), non-biological cover (nobio).

In addition, each of these land surface types is further divided between snow covered and snow-free fractions. For the snow covered fractions, the calculation of the albedo follows a more complex scheme that depends on snow age with a different formulation for the non-biological surface (see details below).

The overall albedo is thus calculated as a surface area-weighted mean of the albedo of each compartment:

Where vegetation albedo ($\alpha^{veg}$) is defined as the combination of bare soil albedo ($alb_\mathsf{bs}$) and leaf albedo for each vegetation type ($alb_\mathsf{leaf}$), weighted by their fractions:

Note that as such the scheme overlooks the effect of vegetation shading the bare soil for sparse canopies, in addition to giving the background for all PFTs the same albedo properties as bare soil fraction of the grid. Snow is allowed to cover both vegetated and non-vegetated areas based on the amount of snow present (see section xxxx). Note also that there is no dependence of the background albedo on soil humidity (i.e, lower albedo when the soil is humid), which may be a limitation in regions with large bare soil fractions.

1.3.2Vegetation albedo¶

We developed a two-stream multilayer matrix-based radiation transfer (RT) solver for the vegetation canopy based on Toon1989. We formulated the RT process as a system of linear equations that considers the interactions of radiative fluxes within and between multiple layers. This approach is particularly well suited for complex vegetation canopies for which their vertical heterogeneity can significantly impact RT. In the two-stream approximation, the basic equations for the upward ($F^\uparrow$) and downward ($F^\downarrow$) diffuse fluxes are given by Meador & Weaver, 1980Toon1989:

The gamma coefficients ($\gamma_1$, $\gamma_2$, $\gamma_3$, and $\gamma_4=1-\gamma_3$) are key parameters in the two-stream approximation that prescribe the scattering and absorption properties of the medium and are thus essential in determining the radiative fluxes within each layer. In particular, $\gamma_1$, and $\gamma_2$ relate to the interactions of upward and downward fluxes, while $\gamma_3$ and $\gamma_4$ represent the source terms associated with scattering and absorption processes. The values of these coefficients can be derived from various approximation methods such as the Eddington approximation, quadrature, and $\delta$-methods Meador & Weaver, 1980. In Eqs. (11), $\tau$ is the optical depth, $\omega$ is the single-scattering albedo, $\pi F_{\odot}$ is the incident flux at the top of the canopy, and $\mu_0$ is the cosine of the solar zenith angle.

When in Eqs. (11) $\gamma_1$ is substituted with $[{1 - \omega \, (1 - \beta)}]/{\bar{\mu}}$, $\gamma_2$ is replaced by ${\omega \, \beta}/{\bar{\mu}}$, $\gamma_3$ is taken as $\beta_0$, $\mu_0$ is redefined as ${1}/{K}$, $\pi F_{\odot}$ is set equal to $K$, and $\tau$ is replaced with $L$ (leaf area index) the two-stream RT equations for the vegetation canopies are obtained Sellers, 1985Dickinson1983Pinty et al., 2006:

Extending the single-layer two-stream model [Eqs. (11)] to multiple layers involves calculating the fluxes within each layer, considering the cumulative optical depths $\tau_c$. In a given layer $l$, the fluxes are expressed as follows:

where

Equations (13) give the analytical solution for Eqs. (11), where $\tau$ is the optical depth within each layer, $\tau_c$ is the cumulative optical depth at the top of each layer (Fig. 3). The coefficients $B_{l1}$ and $B_{l2}$ are derived from the boundary conditions and the continuity at the interfaces between the layers. For numerical stability and to keep all arguments to the exponential terms negative, the variables $A_l$ and $B_l$ as in Eqs. (15) are introduced. This mitigates the risk of numerical overflow or underflow, improving the stability and accuracy of the solver.

In these equations, $\Delta\tau$ represents the total optical depth of the layer (Fig. 3). The downward and upward fluxes then read:

In multilayer matrix-based RT, flux continuity and boundary conditions are maintained using a matrix of coefficients. The continuity equations for upward and downward fluxes at the interfaces of layers $l$ and $l+1$ are given by:

Using Eqs. (16), Eqs. (17) can be expressed as a system of linear equations:

where the coefficients $e_l$, $d_l$, $f_l$, and $g_l$ are defined as:

We derived a tridiagonal system of equations to reduce computational efforts. The following steps outline the derivation process. Multiply the first equation in Eqs. (18) by $d_{l+1}$ and the second by $g_{l+1}$, subtract the resulting equations to eliminate $B_{l+1}$, leading to a new equation relating $A_l$, $B_l$, and $A_{l+1}$. Multiply the second equation in Eqs. (18) by $d_{l+1} e_l - g_{l+1} f_l$ and the first by $f_l$, subtract these results to eliminate $A_{l+1}$, leading to another new equation that relates $A_l$, $B_l$, and $B_{l+1}$. The final system of equations reads:

where the coefficients $\alpha_{l1}$, $\beta_{l1}$, $\gamma_{l1}$ and $\alpha_{l2}$, $\beta_{l2}$, $\gamma_{l2}$ are derived from the original coefficients $e_l$, $d_l$, $f_l$, $g_l$ and the source terms $C_l^{\uparrow}$, $C_l^{\downarrow}$.

The downward diffuse flux at the top of the multilayer structure is equal to any incident diffuse flux. In addition, the upward flux at the bottom of the multilayer structure is equal to the product of the downward flux and the reflectance of the surface ($R_{S}$), i.e. the incident flux reflected at the surface is diffuse. Hence, boundary conditions at the top and bottom boundaries are as follows:

We need to solve a tridiagonal system ($2N-2$ equations) to determine $A_l$, $B_l$, $A_{l+1}$, $B_{l+1}$ [Eqs. (20)] and the related variables such as $F_l^{\downarrow}$ and $F_l^{\uparrow}$. We use the standard method of tridiagonal solving, known as the Thomas algorithm. This algorithm is efficient for solving tridiagonal systems, as it reduces computational complexity compared with general matrix solvers.

Several key parameters need to be precalculated to model RT in the vegetation canopy. These are $G(\mu)$, $K(\mu)$, $\bar{\mu}$, $\omega \, \beta$, and $\omega \, \beta_0$. The function $G(\mu)$, known as the asymmetry factor of the phase function, is given by Myneni et al., 1989:

where $g(\theta)$ and $h(\phi)$ are the probability density functions of leaf normal inclination and azimuth, respectively. Several common functions exist for the probability distribution $g(\theta)\sin \theta$. For instance, the planophile function is given by ${2}/{\pi} \left(1 + \cos2\theta\right)$, while the erectophile function is defined as ${2}/{\pi} \left(1 - \cos2\theta\right)$. The plagiophile function is expressed as ${2}/{\pi} \left(1 - \cos4\theta\right)$, and the extremophile function is represented by ${2}/{\pi} \left(1 + \cos4\theta\right)$. The function is $\sin\theta$ in the spherical case ($G=1/2$) while the uniform function is ${2}/{\pi}$. Once the $G$ function has been determined, the other parameters can be directly derived. These are:

where $r$ is the reflectivity and $t$ is the transmissivity of vegetation (leaf). The last two equations in Eqs. (23) represent the back-scattering parameters ($\beta$, $\beta_0$) for the diffuse and direct beams, respectively, as defined by Pinty et al. (2006). However, these two parameters need to be adjusted for other cases. In Dickinson1983Sellers (1985)Oleson et al. (2010), the back-scattering parameter for the diffuse beam is defined as

where $\Theta$ is the mean leaf inclination angle relative to the horizontal surface. The back-scattering parameter for direct beam, $\beta_0$, reads Dickinson1983Sellers, 1985Oleson et al., 2010:

where $\alpha_{ss}(\mu)$, the single scattering albedo of the canopy, is given by

1.3.3OK: Background and bare soil albedo¶

We combine into a single variable referred as soil albedo i) the background albedo that corresponds to the albedo of the soil under a vegetation cover and that is used as a boundary condition of the vegetation radiative transfer scheme to compute the overall vegetation albedo (see section above) and ii) the bare soil albedo from non-vegetated surfaces (i.e. Bare soil PFT).

Currently, there are two options to define the albedo of the soil. The first option (the default) is a spatially and temporally variable soil albedo, derived from MODIS satellite observations (Ref) using an extensive optimization procedure (see https://orchidas.lsce.ipsl.fr/dev/albedo/optimization.php). A 3-D field (longitude, latitude, time) of soil albedos was derived at a 1° spatial resolution and for a climatology of 12 months. The computation arises from an inverse procedure (similar to Pinty et al. (2011)), where for each 1° grid-cell the monthly soil albedo results from an optimization that accounts for vegetation albedo (see above) and snow albedo (see below). The second option corresponds to the use of a mean value (fixed in time) for each plant functional type that the user can easily prescribe.

1.3.4Snow albedo¶

The snow albedo is computed following the formulation of Chalita & Le Treut (1994) for each PFT, following the equation:

$A_{aged}$ represent albedo of old snow and the sum of $A_{aged}$ and $B_{dec}$ is the albedo of fresh snow. $age_{dec}$ is the time constant of the albedo decay of snow in days and $age_{snow}$ the age of snow parameterized as follows:

where $age_{max}$ is the maximum snow age, $P_{snow}$ is the amount of snowfall during the time interval $dt$ and $\delta_c$ is the solid precipitation constant of snow age. Grid point albedo is weighted according to the fraction of each PFT, distinguishing between forested and grasses and crops PFTs.

For forest PFTs, the albedo is weighted with $f^{veg,max}_{i}$, the maximum cover fraction of a PFT, assuming that even snow-covered vegetation remains visible. For grasses and crops PFTs, the albedo is weighted with $f^{veg}_{i}$, the fraction of vegetation that covers the soil. When $f^{veg}_{i}$ is lower than $f^{veg,max}_{i}$, the albedo of the snow on the baresoil is used. Baresoil snow albedo is weighted by $f_{totbaresoilnotree}$ with

1.3.5Ice sheet and glacier albedo¶

The non-biological compartment currently represents area covered by ice sheet or glaciers. Low surface air temperatures found in this regions slow down the snow metamorphism which is why albedo evolution is specific to these regions. This effect is accounted for with the function $f_{age,nobio}$:

where $\omega1$ and $\omega2$ are tuning constant.

1.4Turbulent transfers¶

1.4.1OK: Main equation for sensible heat, latent heat and momentum fluxes¶

As shown in the surface energy budget in Eq. 6, part of the energy received at the surface is transferred to the atmosphere through turbulent fluxes. These turbulent exchanges are typically represented by three main fluxes: the sensible and latent heat fluxes, which are driven by heat and moisture exchanges between the surface and the atmosphere, and the momentum flux, which represents the drag exerted by the atmosphere on the surface.

Over the grid cell, ORCHIDEE explicitly computes the turbulent fluxes of sensible $F^{H}$ and latent $F^{LE}$ heat (W m$^{-2}$) between the surface and the atmosphere. These fluxes are represented using the diffusive formulation of Budyko (1961) and can be expressed as:

Here, $\rho$ is the air density (kg m$^{-3}$), $c^{cp}$ is the specific heat capacity of air (m$^2$ s$^{-2}$ K$^{-1}$), $\lambda$ is the latent heat of vaporization (or sublimation in the case of a snow-covered surface) (J kg$^{-1}$), $T^{surf}$ is the surface temperature (K) and $q^{sat}(T^{surf})$ is the saturated specific humidity at $T^{surf}$ (kg kg$^{-1}$), $T^{atm}$ and $q^{atm}$ are the atmospheric temperature and specific humidity (K and kg kg$^{-1}$), $\beta$ represents the combined effect of physiological and physical resistances that reduce the actual evapotranspiration relative to the potential evapotranspiration ($F^{ETP}$) (unitless) (see Section 1.6), and $R^{aero}$ is the aerodynamic resistance (s m$^{-1}$). The aerodynamic transport model used to compute the aerodynamic resistance is described in section 1.4.2, while the resistances for sensible and latent heat fluxes are detailed in sections 1.5 and 1.6, respectively.

Although the momentum flux itself is not computed directly in ORCHIDEE, it represents the drag force that generates turbulence and determines the efficiency of heat and moisture transfers through its control on aerodynamic resistances. The zonal ($\tau_x$) and meridional ($\tau_y$) components of the momentum flux can be expressed as:

Where $u^{surf}$ and $v^{surf}$ are the zonal and meridional wind components at the surface (assumed to be zero at the roughness height), $u^{atm}$ and $v^{atm}$ are the corresponding wind components in the atmosphere (m s$^{-1}$), and $R^{aero}$ is the aerodynamic resistance (s m$^{-1}$).

1.4.2OK: Aerodynamic transfer: roughness length and drag coefficients¶

Does this section describe the different options in ORCHIDEE? For example, rough_dyn, use_ratio_z0m_z0h, RATIO_Z0M_Z0H, and use_height_dom. The parameter names should not neceassarily be mentioned but the different options should as they are referred in section 1.

The efficiency of turbulent fluxes introduced in the previous section depends on how the surface interacts with the atmosphere, which is captured through drag coefficients and roughness lengths. When ORCHIDEE is coupled with LMDZ, the drag coefficients are provided by LMDZ, while ORCHIDEE computes the roughness lengths. In offline mode, ORCHIDEE calculates both the drag coefficients and the roughness lengths directly. Although the energy budget is resolved at the grid-cell scale, the drag coefficients are calculated for each PFT, as they depend on vegetation structure, and then aggregated to obtain a grid cell average.

In offline mode, the roughness length is first computed for bare soil $z^{m,bare}_{0}$ (m), as well as the drag coefficients for heat $C^{h,bare}$ and momentum $C^{m,bare}$ (unitless) using the von Karman constant:

Where $f^{snow,veg}$ is the fraction of snow on vegetated area, and $z^{bare}_{0}$ is the bare soil roughness length set to 0.01m.

Where $f^{bare}$ is the total evaporating bare soil fraction, $\kappa=0.41$ is the von Karman constant (unitless), $z^{atm}_{1,max}$ is the maximum height of the first atmospheric layer (set to at least 10m), and $c^{z^{m}_{0},z^{h}_{0}}$ a constant representing the ratio between $z^{m}_{0}$ and $z^{h}_{0}$, set by default to 1 for all PFTs.

Then for each vegetated PFT, several options exist to compute each of these two roughness lengths, which can be combined for a given simulation.

A dynamical computation uses the formulation proposed by Ershadi et al. (2015) based on the initial work of Su et al. (2001) for $z^{m,pft}_{0}$, which is calculated as:

with ${h_{c,pft}}$ the canopy height (m), ${d_{0,pft}}$ the displacement height, assumed equal to $\frac{2}{3}{h_{c,pft}}$, ${k}$ the von Karman’s constant and $\mathsf{\eta}$, the ratio of friction velocity to the wind speed at the top of canopy.

The roughness length for heat is calculated, assuming an exponential relationship between $z^{m,pft}_{0}$ and $z^{h,pft}_{0}$, with the following equation:

with $B$ the Stanton number (see Appendix E of Ershadi et al. (2015) for a full description).

It is also possible to use fixed ratios to define the roughness height for momentum as a fraction of the canopy height, and the roughness height for heat as a fraction of the roughness height for momentum.

Using these values of roughness lengths, the drag coefficients for heat $C^{h,veg}$ and momentum $C^{m,veg}$ (unitless) are computed following:

Where $f^{veg}$ is the fraction of vegetated PFT. Note that in the case of forest PFTs, $f^{veg}$ is set to $f^{veg,max}$ because trees trunks influence the roughness even without leaves, while for grassland PFT, they influence the roughness only during the growing season. $d^{h}$ is the PFT vegetation height, and $c^{z,h}$ is a roughness scaling factor set to 0.0625 for all vegetated PFTs (unitless).

Then, for each none vegetated surface, the drag coefficients for heat $C^{h,nobio}$ and momentum $C^{m,nobio}$ are computed in a similar way:

Where $f^{nobio}$ is the fraction of the non vegetated surface, and $z^{m,nobio}_{0}$ is the roughness height of the non vegetated surface set to 0.001m.

Finally, the grid-cell drag coefficients are obtained by summing the fraction-weighted contributions from bare soil, vegetated PFTs, and non-vegetated surfaces. For vegetated areas, the total is normalized by the effective vegetated fraction to ensure consistent weighting across surface types.

Then, at the grid cell level, the roughness length for heat $z^{h}_{0}$ and for momentum $z^{m}_{0}$ are calculated based on the inversion of the calculation of the drag coefficients aggregated over the grid cell ($C^{h}$ and $C^{m}$):

The zero-plane displacement height $z^{d}$ is also computed as a fraction of average vegetation height:

Where $d^{h,ave}$ is the average vegetation height (m), and $c^{z^{d}}=0.66$ is a factor to calculate the zero-plane displacement height (unitless).

The grid cell effective roughness height $z^{rough}$ is computed as the difference between the average vegetation height and displacement height:

Note that the vertical reference for heights differs between the land and atmospheric components: in ORCHIDEE, heights are expressed relative to the canopy top or to $z^{d}$, while in LMDZ they are referenced from the displacement height. Therefore, $z^{rough}$ provides a consistent means to convert between these reference levels.

The aerodynamic drag coefficient determines the strength of turbulent exchanges between the land surface and the atmosphere. It depends on both surface roughness and atmospheric stability, which is expressed through the Richardson number $Ri$. $Ri$ is defined as the ratio of the two source terms for turbulent kinetic energy: buoyancy and wind shear, which depend on the gradient of temperature near the surface and on the surface wind speed:

Where $z^{atm}_{1}$ is the height of the first atmospheric layer (m), $g$ is the gravitational constant, $u$ is the wind speed ($m\,s^{-1}$), and $T^{v,atm}$ and $T^{v,surf}$ are the virtual temperature of the atmosphere and the surface (K).

Positive $Ri$ values correspond to stable conditions (suppressed turbulence), negative values to unstable conditions (enhanced turbulence), and $Ri=0$ corresponds to neutral conditions.

The neutral drag coefficient $C^{neut}$ (assuming neutral atmospheric stability) is computed as: