The interactions between soil–biosphere–atmosphere land surface model with a multi-energy balance (ISBA-MEB) option in SURFEXv8 – Part 1: Model description

Abstract

Land surface models (LSMs) are pushing towards improved realism owing to an increasing number of observations at the local scale, constantly improving satellite data sets and the associated methodologies to best exploit such data, improved computing resources, and in response to the user community. As a part of the trend in LSM development, there have been ongoing efforts to improve the representation of the land surface processes in the interactions between the soil–biosphere–atmosphere (ISBA) LSM within the EXternalized SURFace (SURFEX) model platform. The force–restore approach in ISBA has been replaced in recent years by multi-layer explicit physically based options for sub-surface heat transfer, soil hydrological processes, and the composite snowpack. The representation of vegetation processes in SURFEX has also become much more sophisticated in recent years, including photosynthesis and respiration and biochemical processes. It became clear that the conceptual limits of the composite soil–vegetation scheme within ISBA had been reached and there was a need to explicitly separate the canopy vegetation from the soil surface. In response to this issue, a collaboration began in 2008 between the high-resolution limited area model (HIRLAM) consortium and Météo-France with the intention to develop an explicit representation of the vegetation in ISBA under the SURFEX platform. A new parameterization has been developed called the ISBA multi-energy balance (MEB) in order to address these issues. ISBA-MEB consists in a fully implicit numerical coupling between a multi-layer physically based snowpack model, a variable-layer soil scheme, an explicit litter layer, a bulk vegetation scheme, and the atmosphere. It also includes a feature that permits a coupling transition of the snowpack from the canopy air to the free atmosphere. It shares many of the routines and physics parameterizations with the standard version of ISBA. This paper is the first of two parts; in part one, the ISBA-MEB model equations, numerical schemes, and theoretical background are presented. In part two (Napoly et al., 2016), which is a separate companion paper, a local scale evaluation of the new scheme is presented along with a detailed description of the new forest litter scheme.

Introduction

Land surface models (LSMs) are based upon fundamental mathematical laws and physics applied within a theoretical framework. Certain processes are modeled explicitly while others use more conceptual approaches. They are designed to work across a large range of spatial scales, so that unresolved scale-dependent processes represented as a function of some grid-averaged state variable using empirical or statistical relationships. LSMs were originally implemented in numerical weather prediction (NWP) and global climate models (GCMs) in order to provide interactive lower boundary conditions for the atmospheric radiation and turbulence parameterization schemes over continental land surfaces. In the past 2 decades, LSMs have evolved considerably to include more biogeochemical and biogeophysical processes in order to meet the growing demands of both the research and the user communities . A growing number of state-of-the-art LSMs, which are used in coupled atmospheric models for operational numerical weather prediction , climate modeling , or both , represent most or all of the following processes: photosynthesis and the associated carbon fluxes, multi-layer soil water and heat transfer, vegetation phenology and dynamics (biomass evolution, net primary production), sub-grid lateral water transfer, river routing, atmosphere–lake exchanges, snowpack dynamics, and near-surface urban meteorology. Some LSMs also include processes describing the nitrogen cycle , groundwater exchanges , aerosol surface emissions , isotopes , and the representation of human impacts on the hydrological cycle in terms of irrigation and ground water extraction , to name a few.

As a part of the trend in LSM development, there have been ongoing efforts to improve the representation of the land surface processes in the Interactions between the soil–biosphere–atmosphere (ISBA) LSM within the EXternalized SURFace (SURFEX; ,) model platform. The original two-layer ISBA force–restore model consists in a single bulk soil layer (generally having a thickness on the order of 50 cm to several meters) coupled to a superficially thin surface composite soil–vegetation–snow layer. Thus, the model simulates fast processes that occur at sub-diurnal timescales, which are pertinent to short-term numerical weather prediction, and it provides a longer-term water storage reservoir, which provides a source for transpiration, a time filter for water reaching a hydro-graphic network, and a certain degree of soil moisture memory in the ground amenable to longer-term forecasts and climate modeling. Additional modifications were made to this scheme over the last decade to include soil freezing , which improved hydrological processes . This scheme was based on the pioneering work of and it has proven its value for coupled land–atmosphere research and applications since its inception. For example, it is currently used for research within the mesoscale non-hydrostatic research model (Meso-NH) . It is also used within the operational high-resolution short-term numerical weather prediction at Météo-France within the limited area model AROME and by HIRLAM countries within the ALADIN–HIRLAM system as the HARMONIE–AROME model configuration . Finally, it is used for climate research within the global climate model (GCM) Action de Researche Petite Echelle Grande Echelle (ARPEGE-climat; , ) and by HIRLAM countries within the ALADIN–HIRLAM system as HARMONIE–AROME and HARMONIE–ALARO Climate configurations .

Rationale for improved vegetation processes

Currently, many LSMs are pushing towards improved realism owing to an increasing number of observations at the local scale, constantly improving satellite data sets and the associated methodologies to best exploit such data, improved computing resources, and in response to the user community via climate services (and seasonal forecasts, drought indexes, etc…). In the SURFEX context, the force–restore approach has been replaced in recent years by multi-layer explicit physically based options for sub-surface heat transfer , soil hydrological processes , and the composite snowpack . These new schemes have recently been implemented in the operational distributed hydro-meteorological hindcast system SAFRAN–ISBA–MODCOU , Meso-NH, and ARPEGE-climat and ALADIN–HIRLAM HARMONIE–AROME and HARMONIE–ALARO Climate configurations. The representation of vegetation processes in SURFEX has also become much more sophisticated in recent years, including photosynthesis and respiration , carbon allocation to biomass pools , and soil carbon cycling . However, for a number of reasons it has also become clear that we have reached the conceptual limits of using of a composite soil–vegetation scheme within ISBA and there is a need to explicitly separate the canopy vegetation from the soil surface:

• in order to distinguish the soil, snow, and vegetation surface temperatures since they can have very different amplitudes and phases in terms of the diurnal cycle, and therefore accounting for this distinction facilitates (at least conceptually) incorporating remote sensing data, such as satellite-based thermal infrared temperatures e.g.,, into such models;

• as it has become evident that the only way to simulate the snowpack beneath forests in a robust and a physically consistent manner (i.e., reducing the dependence of forest snow cover on highly empirical and poorly constrained snow fractional cover parameterizations) and including certain key processes (i.e., canopy interception and unloading of snow) is to include a forest canopy above or buried by the ground-based snowpack e.g.,;

• for accurately modeling canopy radiative transfer, within or below canopy turbulent fluxes and soil heat fluxes;

• to make a more consistent photosynthesis and carbon allocation model (including explicit carbon stores for the vegetation, litter, and soil in a consistent manner);

• to allow the explicit treatment of a ground litter layer, which has a significant impact on ground heat fluxes and soil temperatures (and freezing) and, by extension, the turbulent heat fluxes.

In response to this issue, a collaboration began in 2008 between the high-resolution limited area model (HIRLAM) consortium and Météo-France with the intention to develop an explicit representation of the vegetation in ISBA under the SURFEX platform. A new parameterization has been developed called the ISBA multi-energy balance (MEB) in order to account for all of the above issues.

MEB is based on the classic two-source model for snow-free conditions, which considers explicit energy budgets (for computing fluxes) for the soil and the vegetation, and it has been extended to a three-source model in order to include an explicit representation of snowpack processes and their interactions with the ground and the vegetation. The vegetation canopy is represented using the big-leaf method, which lumps the entire vegetation canopy into a single effective leaf for computing energy budgets and the associated fluxes of heat, moisture, and momentum. One of the first examples of a two-source model designed for atmospheric model studies is , and further refinements to the vegetation canopy processes were added in the years that followed leading to fairly sophisticated schemes, which are similar to those used today e.g.,. The two-source big-leaf approach has been used extensively within coupled regional and global scale land–atmosphere models . In addition, more recently multi-layer vegetation schemes have also been developed for application in GCMs .

ISBA-MEB has been developed taking the same strategy that has been used historically for ISBA: inclusion of the key first-order processes while maintaining a system that has minimal input data requirements and computational cost while being consistent with other aspects of ISBA (with the ultimate goal of being used in coupled operational numerical weather forecast and climate models, and spatially distributed monitoring and hydrological modeling systems). In 2008, one of the HIRLAM partners, the Swedish Meteorological and Hydrological Institute (SMHI), had already developed and applied an explicit representation of the vegetation in the Rossby Centre Regional Climate Model (RCA3) used at SMHI . This representation was introduced into the operational NWP HIRLAMv7.3 system, which became operational in 2010. In parallel, the dynamic vegetation model LJP-GUESS was coupled to RCA3 as RCA–GUESS , making it possible to simulate complex biogeophysical feedback mechanisms in climate scenarios. Since then RCA–GUESS has been applied over Europe , Africa , and the Arctic . The basic principles developed by SMHI have been the foundation since the explicit representation of the vegetation was introduced in ISBA and SURFEX, but now in a more general and consistent way. Implementation of canopy turbulence scheme, longwave radiation transmission function, and snow interception formulations in MEB largely follows the implementation used in RCA3 . In addition, we have taken this opportunity to incorporate several new features into ISBA-MEB compared to the original SMHI scheme:

• a snow fraction that can gradually bury the vegetation vertically thereby transitioning the turbulence coupling from the canopy air space directly to the atmosphere (using a fully implicit numerical scheme);

• the use of the detailed solar radiation transfer scheme that is a multi-layer model that considers two spectral bands, direct and diffuse flux components, and the concept of sunlit and shaded leaves, which was primarily developed to improve the modeling of photosynthesis within ISBA (Carrer et al., 2013);

• a more detailed treatment of canopy snow interception and unloading processes and a coupling with the ISBA physically based multi-layer snow scheme;

• a reformulation of the turbulent exchange coefficients within the canopy air space for stable conditions, such as over a snowpack;

• a fully implicit Jacobean matrix for the longwave fluxes from multiple surfaces (snow, below-canopy snow-free ground surface, vegetation canopy);

• all of the energy budgets are numerically implicitly coupled with each other and with the atmosphere using the coupling method adapted from , which was first proposed by ;

• an explicit forest litter layer model (which also acts as the below-canopy surface energy budget when litter covers the soil).

This paper is the first of two parts: in part one, the ISBA-MEB model equations, numerical schemes, and theoretical background are presented. In part two, a local-scale evaluation of the new scheme is presented along with a detailed description of the new forest litter scheme . An overview of the model is given in the next section, followed by conclusions.

Model description

Description of the patches for the natural land surface sub-grid tile. The values for the 19-class option are shown in the leftmost three columns, and those for the 12-class option are shown in the rightmost three columns (the name and description are only given if they differ from the 19-class values). MEB can currently be activated for the forest classes: 4–6 (for both the 12- and 19-class options), and 13–17.

Index	Name	Description	Index	Name	Description
1	NO	Bare soil	1
2	ROCK	Rock	2
3	SNOW	Permanent snow or ice	3
4	TEBD	Temperate broad leaf	4	TREE	Broad leaf
5	BONE	Boreal evergreen needle leaf	5	CONI	Evergreen needle leaf
6	TRBE	Tropical evergreen broad leaf	6	EVER	Evergreen broad leaf
7	C3	C3 crops	7
8	C4	C4 crops	8
9	IRR	Irrigated crops	9
10	GRAS	Temperate grassland	10
11	TROG	Tropical grassland	11
12	PARK	Bog, park, garden	12
13	TRBD	Tropical broad leaf
14	TEBE	Temperate evergreen broad leaf
15	TENE	Temperate evergreen needle leaf
16	BOBD	Boreal broad leaf
17	BOND	Boreal needle leaf
18	BOGR	Boreal grassland
19	SHRB	Shrubs

SURFEX uses the tile approach for the surface, and separate physics modules are used to compute surface–atmosphere exchange for oceans or seas, lakes, urbanized areas, and the natural land surface . The ISBA LSM is used for the latter tile, and the land surface is further split into upwards of 12 or 19 patches (refer to Table ), which represent the various land cover and plant functional types. Currently, forests make up eight patches for the 19-class option, and three for the 12-class option. The ISBA-MEB (referred to hereafter simply as MEB) option can be activated for any number of the forest patches. By default, MEB is coupled to the multi-layer soil (ISBA-DF: explicit DiFfusion equation for heat and Richard's equation for soil water flow; ,; , ) and snow (ISBA-ES: multi-layer Explicit Snow processes with 12 layers by default; , , , ) schemes. These schemes have been recently updated to include improved physics and increased layering (14 soil layers by default). MEB can also be coupled to the simple three-layer soil force–restore (3-L) option in order to be compatible with certain applications, which have historically used 3-L, but by default, it is coupled with ISBA-DF since the objective is to move towards a less conceptual LSM.

A schematic representation of the turbulent aerodynamic resistance, $R_{\mathrm{a}}$, pathways for ISBA-MEB. The prognostic temperature, liquid water, and liquid water equivalent variables are shown. The canopy air diagnostic variables are enclosed by the red-dashed circle. The ground-based snowpack is indicated using turquoise, the vegetation canopy is shaded green, and ground layers are colored dark brown at the surface, fading to white with depth. Atmospheric variables (lowest atmospheric model or observed reference level) are indicated using the a subscript. The ground snow fraction, $p_{\mathrm{ng}}$, and canopy-snow-cover fraction, $p_{\mathrm{n}\mathit{\alpha }}$, are indicated.

[Figure omitted. See PDF]

A schematic diagram illustrating the various resistance pathways corresponding to the turbulent fluxes for the three fully (implicitly) coupled surface energy budgets is shown in Fig. . The water budget prognostic variables are also indicated. Note that the subscripts, which are used to represent the different prognostic and diagnostic variables and the aerodynamic resistance pathways, are summarized in Table . The canopy bulk vegetation layer is represented using green, the canopy-intercepted snow and ground-based snowpack are shaded using turquoise, and the ground layers are indicated using dark brown at the surface, which fade to white with increasing depth.

There are six aerodynamic resistance, $R_{\mathrm{a}}$ (s${}^{-\mathrm{1}}$), pathways defined as being between (i) the non-snow-buried vegetation canopy and the canopy air, $R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}$, (ii) the non-snow-buried ground surface (soil or litter) and the canopy air, $R_{\mathrm{a}\mathrm{g}-\mathrm{c}}$, (iii) the snow surface and the canopy air, $R_{\mathrm{a}\mathrm{n}-\mathrm{c}}$, (iv) the ground-based snow-covered part of the canopy and the canopy air, $R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}$, (v) the canopy air with the overlying atmosphere, $R_{\mathrm{a}\mathrm{c}-\mathrm{a}}$), and (vi) the ground-based snow surface (directly) with the overlying atmosphere, $R_{\mathrm{a}\mathrm{n}-\mathrm{a}}$. Previous papers describing ISBA expressed heat fluxes using a dimensionless heat and mass exchange coefficient, $C_{\mathrm{H}}$; however, for the new MEB option, it is more convenient to express the different fluxes using resistances (s m${}^{-\mathrm{1}}$), which are related to the exchange coefficient as $R_{\mathrm{a}}=\mathrm{1}/\left(V_{\mathrm{a}},,C_{\mathrm{H}}\right)$, where $V_{\mathrm{a}}$ represents the wind speed at the atmospheric forcing level (indicated by using the subscript a) in m s${}^{-\mathrm{1}}$.

The surface energy budgets are formulated in terms of prognostic equations governing the evolutions of the bulk vegetation canopy, $T_{\mathrm{v}}$, the snow-free ground surface (soil or litter), $T_{\mathrm{g}}$, and the ground-based snowpack, $T_{\mathrm{n}}$ (K). The prognostic hydrological variables consist of the liquid soil water content, $W_{\mathrm{g}}$, equivalent water content of ice, $W_{\mathrm{gf}}$, snow water equivalent (SWE), $W_{\mathrm{n}}$, vegetation canopy-intercepted liquid water, $W_{\mathrm{r}}$, and intercepted snow, $W_{\mathrm{rn}}$ (kg m${}^{-\mathrm{2}}$). The diagnosed canopy air variables, which are determined implicitly during the simultaneous solution of the energy budgets, are enclosed within the red-dashed circle and represent the canopy air specific humidity, $q_{\mathrm{c}}$ (kg kg${}^{-\mathrm{1}}$), air temperature $T_{\mathrm{c}}$, and wind speed $V_{\mathrm{c}}$. The ground surface specific humidity is represented by $q_{\mathrm{g}}$. The surface snow cover fraction area is represented by $p_{\mathrm{ng}}$ while the fraction of the canopy buried by the ground-based snowpack is defined as $p_{\mathit{\alpha }\mathrm{n}}$. The snowpack has $N_{\mathrm{n}}$ layers, while the number of soil layers is defined as $N_{\mathrm{g}}$ where $k$ is the vertical index (increasing from 1 at the surface downward). The ground and snowpack uppermost layer temperatures correspond to those used for the surface energy budget (i.e., $k=\mathrm{1}$).

Snow fractions

Subscripts used to represent the prognostic and diagnostic variables. In addition, the symbols used to represent the aerodynamic resistance pathways (between the two elements separated by the dash) are also shown (refer also to Fig. ). These symbols are used throughout the text.

Subscript	Name	Description
v	Vegetation	Bulk canopy layer
g	Ground	Temperature or liquid water (for $N_{\mathrm{g}}$ layers)
gf	Ground	Frozen water (for $N_{\mathrm{g}}$ layers)
a	Atmosphere	At the lowest atmospheric or forcing level
c	Canopy air space	Diagnosed variables
n	Ground-based snowpack	For $N_{\mathrm{n}}$ layers
ng	Ground-based snowpack	Fractional ground snow coverage
$\mathit{\alpha }$ n	Ground-based snowpack	Fractional vegetation snow coverage
r	Interception reservoir	Intercepted rain and snow meltwater
rn	Interception reservoir	Intercepted snow and frozen meltwater or rain
vg-c	Aerodynamic resistance	Non-snow-buried vegetation canopy and canopy air
g-c	Aerodynamic resistance	Non-snow-buried ground surface and canopy air
n-c	Aerodynamic resistance	Snow surface and canopy air
vn-c	Aerodynamic resistance	Ground-based snow-covered canopy and canopy air
c-a	Aerodynamic resistance	Canopy air with overlying atmosphere
n-a	Aerodynamic resistance	Ground-based snow surface and overlying atmosphere

Snow is known to have a significant impact on heat conduction fluxes, owing to its relatively high insulating properties. In addition, it can significantly reduce turbulent transfer owing to reduced surface roughness, and it has a relatively large surface albedo thereby impacting the surface net radiation budget. Thus, the parameterization of its areal coverage turns out to be a critical aspect of LSM modeling of snowpack–atmosphere interactions and sub-surface soil and hydrological processes. The fractional ground coverage by the snowpack is defined as $$p_{\mathrm{ng}}=W_{\mathrm{n}}/W_{\mathrm{n},\mathrm{crit}}\left(\mathrm{0},\le ,p_{\mathrm{ng}},\le ,\mathrm{1}\right),$$ where currently the default value is $W_{\mathrm{n},\mathrm{crit}}=\mathrm{1}$ (kg m${}^{-\mathrm{2}}$). Note that this is considerably lower than the previous value of 10 kg m${}^{-\mathrm{2}}$ used in ISBA , but this value has been shown to improve the ground soil temperatures, using an explicit snow scheme within ISBA .

The fraction of the vegetation canopy, which is buried by ground-based snow, is defined as $$p_{\mathrm{n}\mathit{\alpha }}=\left(D_{\mathrm{n}},-,z_{\mathrm{hv},\mathrm{b}}\right)/\left(z_{\mathrm{hv}},-,z_{\mathrm{hv},\mathrm{b}}\right)\left(\mathrm{0},\le ,p_{\mathrm{n}\mathit{\alpha }},\le ,\mathrm{1}\right),$$ where $D_{\mathrm{n}}$ is the total ground-based snowpack depth (m) and $z_{\mathrm{hvb}}$ represents the base of the vegetation canopy (m) (see Fig. ), which is currently defined as $$z_{\mathrm{hvb}}=a_{\mathrm{hv}}\left(z_{\mathrm{hv}},-,z_{\mathrm{hv},\min }\right)\left(z_{\mathrm{hvb}},\ge ,\mathrm{0}\right),$$ where $a_{\mathrm{hv}}=0.2$ and the effective canopy base height is set to $z_{\mathrm{hv},\min }=\mathrm{2}$ (m) for forests. The foliage distribution should be reconsidered in further development since literature suggests, e.g., , that the foliage is not symmetrically distributed in the crown but skewed upward.

Energy budget

The coupled energy budget equations for a three-source model can be expressed for a single bulk canopy, a ground-based snowpack, and a underlying ground surface as $$\begin{array}{cc}& {\displaystyle }{\mathcal{C}}_{\mathrm{v}}{\displaystyle \frac{\partial T_{\mathrm{v}}}{\partial t}}=& {\displaystyle }{R}_{\mathrm{n}\mathrm{v}}-H_{\mathrm{v}}-L{E}_{\mathrm{v}}+L_{\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{v}},{\displaystyle }{\mathcal{C}}_{\mathrm{g},\mathrm{1}}{\displaystyle \frac{\partial T_{\mathrm{g},\mathrm{1}}}{\partial t}}=& {\displaystyle }\left(\mathrm{1},-,p_{\mathrm{ng}}\right)\left(R_{\mathrm{n}\mathrm{g}},-,H_{\mathrm{g}},-,{\mathrm{LE}}_{\mathrm{g}}\right)\\ {\displaystyle }& {\displaystyle }+p_{\mathrm{ng}}\left(G_{\mathrm{gn}},+,{\mathit{\tau }}_{\mathrm{n},N_{\mathrm{n}}},{\mathrm{SW}}_{\mathrm{net},\mathrm{n}}\right)\\ & {\displaystyle }& {\displaystyle }-G_{\mathrm{g},\mathrm{1}}+L_{\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{g},\mathrm{1}},{\displaystyle }{\mathcal{C}}_{\mathrm{n},\mathrm{1}}{\displaystyle \frac{\partial T_{\mathrm{n},\mathrm{1}}}{\partial t}}=& {\displaystyle }{R}_{\mathrm{n}\mathrm{n}}-H_{\mathrm{n}}-{\mathrm{LE}}_{\mathrm{n}}-{\mathit{\tau }}_{\mathrm{n},\mathrm{1}}{\mathrm{SW}}_{\mathrm{net},\mathrm{n}}+{\mathit{\xi }}_{\mathrm{n},\mathrm{1}}\\ & {\displaystyle }& {\displaystyle }-G_{\mathrm{n},\mathrm{1}}+L_{\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{n},\mathrm{1}},\end{array}$$ where $T_{\mathrm{g},\mathrm{1}}$ is the uppermost ground (surface soil or litter layer) temperature, $T_{\mathrm{n},\mathrm{1}}$ is the surface snow temperature, and $T_{\mathrm{v}}$ is the bulk canopy temperature (K). Note that the subscript 1 indicates the uppermost layer or the base of the layer (for fluxes) for the soil and snowpack. All of the following flux terms are expressed in W m${}^{-\mathrm{2}}$. The sensible heat fluxes are defined between the canopy air space and the vegetation $H_{\mathrm{v}}$, the snow-free ground $H_{\mathrm{g}}$, and the ground-based snowpack $H_{\mathrm{n}}$. In an analogous fashion to the sensible heat flux, the latent heat fluxes are defined for the vegetation canopy $E_{\mathrm{v}}$, the snow-free ground $E_{\mathrm{g}}$, and the ground-based snowpack $E_{\mathrm{n}}$. The net radiation fluxes are defined for the vegetation canopy, ground, and snowpack as $R_{\mathrm{n}\mathrm{v}}$, $R_{\mathrm{n}\mathrm{g}}$, and $R_{\mathrm{n}\mathrm{n}}$, respectively. Note that part of the incoming shortwave radiation is transmitted through the uppermost snow layer, and this energy loss is expressed as ${\mathit{\tau }}_{\mathrm{n},\mathrm{1}}{\mathrm{SW}}_{\mathrm{net},\mathrm{n}}$, where $\mathit{\tau }$ is the dimensionless transmission coefficient. The conduction fluxes between the uppermost ground layer and the underlying soil and the analogue for the snowpack are defined as $G_{\mathrm{g},\mathrm{1}}$ and $G_{\mathrm{n},\mathrm{1}}$, respectively. The conduction flux between the base of the snowpack and the ground surface is defined as $G_{\mathrm{gn}}$. The last term on the right-hand side (RHS) of Eq. (), ${\mathit{\xi }}_{\mathrm{n},\mathrm{1}}$, represents the effective heating or cooling of a snowpack layer caused by exchanges in enthalpy between the surface and sub-surface model layers when the vertical grid is reset (the snow model grid-layer thicknesses vary in time).

The ground-based snow fraction is defined as $p_{\mathrm{ng}}$. Note that certain terms of Eq. () are multiplied by $p_{\mathrm{ng}}$ to make them patch relative (or grid box relative in the case of single-patch mode) since the snow can potentially cover only part of the patch. Within the snow module itself, the notion of $p_{\mathrm{ng}}$ is not used (the computations are snow relative). But note that when simultaneously solving the coupled equations Eqs. ()–(), Eq. () must be multiplied by $p_{\mathrm{ng}}$ since again, snow only covers a fraction of the area: further details are given in Appendices  and . The formulation for $p_{\mathrm{ng}}$ is described in Sect. .

A schematic sketch illustrating the role of $p_{\mathrm{n}\mathit{\alpha }}$, the fraction of the vegetation canopy, which is buried by ground-based snow. In panel (a), the snow is well below the canopy base, $z_{\mathrm{hvb}}$, resulting in $p_{\mathrm{n}\mathit{\alpha }}=\mathrm{0}$, and the snow has no direct energy exchange with the atmosphere. In panel (b), the canopy is partly buried by snow ($\mathrm{0}<p_{\mathrm{n}\mathit{\alpha }}<\mathrm{1}$) and the snow has energy exchanges with both the canopy air and the atmosphere. In panel (c), the canopy is fully buried by snow ($p_{\mathrm{n}\mathit{\alpha }}=\mathrm{1}$) and the snow has energy exchange only with the atmosphere, whereas the soil and canopy only exchange with the canopy air space ($p_{\mathrm{ng}}<\mathrm{1}$). Finally, in panel (d), both $p_{\mathrm{ng}}=\mathrm{1}$ and $p_{\mathrm{n}\mathit{\alpha }}=\mathrm{1}$ so that the only exchanges are between the snow and the atmosphere.

[Figure omitted. See PDF]

The phase change terms (freezing less melting: expressed in kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) for the snow water equivalent intercepted by the vegetation canopy, the uppermost ground layer, and the uppermost snowpack layer are represented by ${\mathrm{\Phi }}_{\mathrm{v}}$, ${\mathrm{\Phi }}_{\mathrm{g},\mathrm{1}}$, and ${\mathrm{\Phi }}_{\mathrm{n},\mathrm{1}}$, respectively, and $L_{\mathrm{f}}$ represents the latent heat of fusion (J kg${}^{-\mathrm{1}}$). The computation of ${\mathrm{\Phi }}_{\mathrm{g},\mathrm{1}}$ uses the Gibbs free-energy method , ${\mathrm{\Phi }}_{\mathrm{n},\mathrm{1}}$ is based on available liquid for freezing or cold content for freezing , and ${\mathrm{\Phi }}_{\mathrm{v}}$ is described herein (see Eq. ). Note that all of the phase change terms are computed as adjustments to the surface temperatures (after the fluxes have been computed); therefore, only the energy storage terms are modified directly by phase changes for each model time step.

The surface ground, snow, and vegetation effective heat capacities, ${\mathcal{C}}_{\mathrm{g},\mathrm{1}}$, ${\mathcal{C}}_{\mathrm{v}}$, and ${\mathcal{C}}_{\mathrm{n},\mathrm{1}}$ (J m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$) are defined, respectively, as $$\begin{array}{}& {\displaystyle }& {\displaystyle }{\mathcal{C}}_{\mathrm{g},\mathrm{1}}=\mathrm{\Delta }{z}_{\mathrm{g},\mathrm{1}}{c}_{\mathrm{g},\mathrm{1}}& {\displaystyle }& {\displaystyle }{\mathcal{C}}_{\mathrm{v}}=C_{\mathrm{vb}}+C_{\mathrm{i}}{W}_{\mathrm{r},\mathrm{n}}+C_{\mathrm{w}}{W}_{\mathrm{r}},& {\displaystyle }& {\displaystyle }{\mathcal{C}}_{\mathrm{n},\mathrm{1}}=D_{\mathrm{n},\mathrm{1}}{c}_{\mathrm{n},\mathrm{1}},\end{array}$$ where $C_{\mathrm{i}}$ and $C_{\mathrm{w}}$ are the specific heat capacities for solid ($2.106\times {10}^{\mathrm{3}}$ J kg${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}}$) and liquid water ($4.218\times {10}^{\mathrm{3}}$ J kg${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}}$), respectively. The uppermost ground-layer thickness is $\mathrm{\Delta }{z}_{\mathrm{g},\mathrm{1}}$ (m), and the corresponding heat capacity of this layer is defined as $c_{\mathrm{g}\mathrm{\hspace{0.17em}1}}$ (J m${}^{-\mathrm{3}}$ K${}^{-\mathrm{1}}$). The uppermost soil layer ranges between 0.01 and 0.03 m for most applications so that the interactions between surface fluxes and fast temperature changes in the surface soil layer can be represented. There are two options for modeling the thermal properties of the uppermost ground layer. First, they can be defined using the default ISBA configuration for a soil layer with parameters based on soil texture properties, which can also incorporate the thermal effects of soil organics . The second option, which is the default when using MEB, is to model the uppermost ground layer as forest litter. The ground surface in forest regions is generally covered by a litter layer consisting of dead leaves and or needles, branches, fruit, and other organic material. Some LSMs have introduced parameterizations for litter , but the approach can be very different from one case to another depending on their complexity. The main goal of this parameterization within MEB is to account for the generally accepted first-order energetic and hydrological effects of litter; this layer is generally accepted to have a strong insulating effect owing to its particular thermal properties (leading to a relatively low thermal diffusivity), it causes a significant reduction of ground evaporation (capillary rise into this layer is negligible), and it constitutes an interception reservoir for liquid water, which can also lose water by evaporation. See for a detailed description of this scheme and its impact on the surface energy budget.

The canopy is characterized by low heat capacity, which means that its temperature responds fast to changes in fluxes. Thus, to realistically simulate diurnal variations in 2 m temperature this effect must be accounted for. defined the value as being the heat capacity of 0.2 kg m${}^{-\mathrm{2}}$ of water per unit leaf area index (m${}^{\mathrm{2}}$ m${}^{-\mathrm{2}}$). This results in values on the order of $\mathrm{1}\times {10}^{\mathrm{4}}$ J m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$ for forest canopies in general. For local-scale simulations, $C_{\mathrm{vb}}$ can be defined based on observational data. In spatially distributed simulations (or when observational data is insufficient), $C_{\mathrm{vb}}=0.2/C_{\mathrm{V}}$ where the vegetation thermal inertia $C_{\mathrm{V}}$ is defined as a function of vegetation class by the SURFEX default physiographic database ECOCLIMAP . Note that $C_{\mathrm{V}}$ has been determined for the composite soil–vegetation scheme, and the factor 0.2 is used to reduce this value to be more representative of vegetation and on the order of the value discussed by . Numerical tests have shown that using this value, the canopy heat storage is on the order of 10 W m${}^{-\mathrm{2}}$ at mid-day for a typical mid-latitude summer day for a forest. The minimum vegetation heat capacity value is limited at $\mathrm{1}\times {10}^{\mathrm{4}}$ (J m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$) in order to model, in a rather simple fashion, the thermal inertia of stems, branches, trunks, etc. The contributions from intercepted snow and rain are incorporated, where $W_{\mathrm{r},\mathrm{n}}$ and $W_{\mathrm{r}}$ (kg m${}^{-\mathrm{2}}$) represent the equivalent liquid water content of intercepted canopy snow and liquid water, respectively.

The uppermost snow-layer thickness is $D_{\mathrm{n},\mathrm{1}}$ (m), and the corresponding heat capacity is represented by $c_{\mathrm{n},\mathrm{1}}$ . Note that $D_{\mathrm{n},\mathrm{1}}$ is limited to values no larger than several centimeters in order to model a reasonable thermal inertia (i.e., in order to represent the diurnal cycle) in a fashion analogous to the soil. For more details, see .

The numerical solution of the surface energy budget, sub-surface soil and snow temperatures, and the implicit numerical coupling with the atmosphere is described in Appendix .

Turbulent fluxes

In this section, the turbulent heat and water vapor fluxes in Eqs. ()–() are described.

Sensible heat fluxes

The MEB sensible heat fluxes are defined as $$\begin{array}{}& {\displaystyle }{H}_{\mathrm{v}}=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}{\displaystyle \frac{\left({\mathcal{T}}_{\mathrm{v}},-,{\mathcal{T}}_{\mathrm{c}}\right)}{R_{\mathrm{a}\mathrm{v}-\mathrm{c}}}},& {\displaystyle }{H}_{\mathrm{g}}=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}{\displaystyle \frac{\left({\mathcal{T}}_{\mathrm{g}},-,{\mathcal{T}}_{\mathrm{c}}\right)}{R_{\mathrm{a}\mathrm{g}-\mathrm{c}}}},& {\displaystyle }{H}_{\mathrm{n}}=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}\left[\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right),,{\displaystyle \frac{\left({\mathcal{T}}_{\mathrm{n}},-,{\mathcal{T}}_{\mathrm{c}}\right)}{R_{\mathrm{a}\mathrm{n}-\mathrm{c}}}},,+,,p_{\mathrm{n}\mathit{\alpha }},,{\displaystyle \frac{\left({\mathcal{T}}_{\mathrm{n}},-,{\mathcal{T}}_{\mathrm{a}}\right)}{R_{\mathrm{a}\mathrm{n}-\mathrm{a}}}}\right],& {\displaystyle }{H}_{\mathrm{c}}=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}{\displaystyle \frac{\left({\mathcal{T}}_{\mathrm{c}},-,{\mathcal{T}}_{\mathrm{a}}\right)}{R_{\mathrm{a}\mathrm{c}-\mathrm{a}}}},& {\displaystyle }H=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}\left[\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }},,p_{\mathrm{ng}}\right),,{\displaystyle \frac{\left({\mathcal{T}}_{\mathrm{c}},-,{\mathcal{T}}_{\mathrm{a}}\right)}{R_{\mathrm{a}\mathrm{c}-\mathrm{a}}}},,+,,p_{\mathrm{n}\mathit{\alpha }},,p_{\mathrm{ng}},,{\displaystyle \frac{\left({\mathcal{T}}_{\mathrm{n}},-,{\mathcal{T}}_{\mathrm{a}}\right)}{R_{\mathrm{a}\mathrm{n}-\mathrm{a}}}}\right],\end{array}$$ where ${\mathit{\rho }}_{\mathrm{a}}$ represents the lowest atmospheric layer average air density (kg m${}^{-\mathrm{3}}$). The sensible heat fluxes appear in the surface energy budget equations (Eqs. –). The sensible heat flux from the ground-based snowpack (Eq. ) is partitioned by the fraction of the vegetation, which is buried by the ground-based snowpack $p_{\mathrm{n}\mathit{\alpha }}$, between an exchange between the canopy air space, and the overlying atmosphere (Eq. ). The heat flux between the overlaying atmosphere and the canopy air space is represented by $H_{\mathrm{c}}$, and it is equivalent to the sum of the fluxes between the different energy budgets and the canopy air space. The total flux exchange between the overlying atmosphere and the surface (as seen by the atmosphere) is defined by $H$. It is comprised of two components: the heat exchange between the overlying atmosphere and the canopy air space and the part of the ground-based snowpack that is burying the vegetation. This method has been developed to model the covering of low vegetation canopies by a ground-based snowpack. Finally, the final fluxes for the given patch are aggregated using $p_{\mathrm{ng}}$ and $p_{\mathrm{n}\mathit{\alpha }}$: the full expressions are given in Appendix .

The thermodynamic variable ($\mathcal{T}$: J kg${}^{-\mathrm{1}}$) is linearly related to temperature as $${\mathcal{T}}_x={\mathcal{B}}_x+{\mathcal{A}}_x{T}_x,$$ where $x$ corresponds to one of the three surface temperatures ($T_{\mathrm{g}}$, $T_{\mathrm{v}}$, or $T_{\mathrm{n}}$), canopy air temperature, $T_{\mathrm{c}}$, or the overlying atmospheric temperature, $T_{\mathrm{a}}$. The definitions of ${\mathcal{A}}_x$ and ${\mathcal{B}}_x$ depend on the atmospheric variable in the turbulent diffusion scheme and are usually defined to cast $\mathcal{T}$ in the form of dry static energy, or potential temperature and are determined by the atmospheric model in coupled mode (see Appendix ). The total canopy aerodynamic resistance is comprised of snow buried, $R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}$, and non-snow buried, $R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}$, resistances from $$R_{\mathrm{a}\mathrm{v}-\mathrm{c}}={\left[{\displaystyle \frac{\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right)p_{\mathrm{ng}}}{R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}}},,+,,{\displaystyle \frac{\left(\mathrm{1},-,p_{\mathrm{ng}}\right)}{R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}}}\right]}^{-\mathrm{1}}.$$ The separation of the resistances is done to mainly account for differences in the roughness length between the buried and non-covered parts of the vegetation canopy; therefore, the primary effect of snow cover is to increase the resistance relative to a snow-free surface assuming the same temperature gradient owing to a lower surface roughness, and thus $R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}\ge R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}$. The formulation also provides a continuous transition to the case of vanishing canopy turbulent fluxes as the canopy becomes entirely buried (as $p_{\mathrm{n}\mathit{\alpha }}\to \mathrm{1}$). In this case, the energy budget equations collapse into a simple coupling between the snow surface and the overlying atmosphere, and the ground energy budget simply consists in heat conduction between the ground surface and the snowpack base. The formulations of the resistances between the different surfaces and the canopy airspace and the overlying atmosphere are described in detail in Sect. . The canopy air temperature, which is needed by different physics routines, is diagnosed by combining Eqs. ()–() and solving for ${\mathcal{T}}_{\mathrm{c}}$ and using Eq. () to determine $T_{\mathrm{c}}$ (see Appendix  for details).

Water vapor fluxes

The MEB water vapor fluxes are expressed as $$\begin{array}{cc}& {\displaystyle }{E}_{\mathrm{v}}=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}{h}_{\mathrm{sv}}{\displaystyle \frac{\left(q_{\mathrm{sat}\mathrm{v}},-,q_{\mathrm{c}}\right)}{R_{\mathrm{a}\mathrm{v}-\mathrm{c}}}},& {\displaystyle }{E}_{\mathrm{g}}=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}{\displaystyle \frac{\left(q_{\mathrm{g}},-,q_{\mathrm{c}}\right)}{R_{\mathrm{a}\mathrm{g}-\mathrm{c}}}},& {\displaystyle }{E}_{\mathrm{n}}=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}{h}_{\mathrm{sn}}\left[\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right),,{\displaystyle \frac{\left(q_{\mathrm{sati}\mathrm{n}},-,q_{\mathrm{c}}\right)}{R_{\mathrm{a}\mathrm{n}-\mathrm{c}}}},,+,,p_{\mathrm{n}\mathit{\alpha }},,{\displaystyle \frac{\left(q_{\mathrm{sati}\mathrm{n}},-,q_{\mathrm{a}}\right)}{R_{\mathrm{a}\mathrm{n}-\mathrm{a}}}}\right],& {\displaystyle }{E}_{\mathrm{c}}=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}{\displaystyle \frac{\left(q_{\mathrm{c}},-,q_{\mathrm{a}}\right)}{R_{\mathrm{a}\mathrm{c}-\mathrm{a}}}},{\displaystyle }E=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}\left[\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }},,p_{\mathrm{ng}}\right),,{\displaystyle \frac{\left(q_{\mathrm{c}},-,q_{\mathrm{a}}\right)}{R_{\mathrm{a}\mathrm{c}-\mathrm{a}}}}\right.\\ & {\displaystyle }& {\displaystyle }\left..,+,,p_{\mathrm{n}\mathit{\alpha }},,p_{\mathrm{ng}},,h_{\mathrm{sn}},,{\displaystyle \frac{\left(q_{\mathrm{sati}\mathrm{n}},-,q_{\mathrm{a}}\right)}{R_{\mathrm{a}\mathrm{n}-\mathrm{a}}}}\right].\end{array}$$ The vapor flux between the canopy air and the overlying atmosphere is represented by $E_{\mathrm{c}}$, and the total vapor flux exchanged with the overlying atmosphere is defined as $E$. The specific humidity (kg kg${}^{-\mathrm{1}}$) of the overlying atmosphere is represented by $q_{\mathrm{a}}$, whereas $q_{\mathrm{sat}}$ and $q_{\mathrm{sati}}$ represent the specific humidity at saturation over liquid water and ice, respectively. For the surface specific humidities at saturation, the convention $q_{\mathrm{sat}x}=q_{\mathrm{sat}}\left(T_x\right)$ is used. The same holds true for saturation over ice so that $q_{\mathrm{sati}\mathrm{n}}=q_{\mathrm{sati}}\left(T_{\mathrm{n}}\right)$. The canopy air specific humidity $q_{\mathrm{c}}$ is diagnosed assuming that $E_{\mathrm{c}}$ is balanced by the vapor fluxes between the canopy air and each of the three surfaces considered (the methodology for diagnosing the canopy air thermal properties is described in Appendix , Sect. ). The effective ground specific humidity is defined as $$q_{\mathrm{g}}=h_{\mathrm{sg}}{q}_{\mathrm{sat}\mathrm{g}}+\left(\mathrm{1},+,h_{\mathrm{a}}\right)q_{\mathrm{c}},$$ where the humidity factors are defined as $$\begin{array}{}& {\displaystyle }{h}_{\mathrm{sg}}=& {\displaystyle }{\mathit{\delta }}_{\mathrm{g}}{h}_{\mathrm{ug}}\left(\mathrm{1},-,p_{\mathrm{gf}}\right)\left({\displaystyle \frac{L_{\mathrm{v}}}{L}}\right)+{\mathit{\delta }}_{\mathrm{gf}}{h}_{\mathrm{ugf}}{p}_{\mathrm{gf}}\left({\displaystyle \frac{L_{\mathrm{s}}}{L}}\right),& {\displaystyle }{h}_{\mathrm{a}}=& {\displaystyle }{\mathit{\delta }}_{\mathrm{g}}\left(\mathrm{1},-,p_{\mathrm{gf}}\right)\left({\displaystyle \frac{L_{\mathrm{v}}}{L}}\right)+{\mathit{\delta }}_{\mathrm{gf}}{p}_{\mathrm{gf}}\left({\displaystyle \frac{L_{\mathrm{s}}}{L}}\right).\end{array}$$ The latent heats of fusion and vaporization are defined as $L_{\mathrm{s}}$ and $L_{\mathrm{v}}$ (J kg${}^{-\mathrm{1}}$), respectively. The fraction of the surface layer that is frozen, $p_{\mathrm{gf}}$, is simply defined as the ratio of the liquid water equivalent ice content to the total water content. The average latent heat $L$ is essentially a normalization factor that ranges between $L_{\mathrm{s}}$ and $L_{\mathrm{v}}$ as a function of snow cover and surface soil ice (see Appendix ). The soil coefficient ${\mathit{\delta }}_{\mathrm{g}}$ in Eqs. ()–() is defined as $${\mathit{\delta }}_{\mathrm{g}}=\left({\displaystyle \frac{R_{\mathrm{a}\mathrm{g}-\mathrm{c}}}{R_{\mathrm{a}\mathrm{g}-\mathrm{c}}+R_{\mathrm{g}}}}\right){\mathit{\delta }}_{\mathrm{gcor}},$$ where the soil resistance $R_{\mathrm{g}}$ is defined by Eq. (). Note that the composite version of ISBA did not include an explicit soil resistance term, and therefore this also represents a new addition to the model. This term was found to further improve results for bare-soil evaporation within MEB, and its inclusion is consistent with other similar multi-source models e.g.,. See Sect.  for further details. The delta function ${\mathit{\delta }}_{\mathrm{gcor}}$ is a numerical correction term that is required owing to the linearization of $q_{\mathrm{sat}\mathrm{g}}$ and is unity unless both $h_{\mathrm{ug}}{q}_{\mathrm{sat}\mathrm{g}}<q_{\mathrm{c}}$ and $q_{\mathrm{sat}\mathrm{g}}>q_{\mathrm{c}}$, in which case it is set to zero. The surface ground humidity factor is defined using the standard ISBA formulation from as $$h_{\mathrm{ug}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{1},,-,,\text{cos},\left({\displaystyle \frac{w_{\mathrm{g},\mathrm{1}}}{w_{\mathrm{fc},\mathrm{1}}^{\ast }}},\mathit{\pi }\right)\right]\left(\mathrm{0},\le ,h_{\mathrm{ug}},\le ,\mathrm{1}\right).$$ In the case of condensation ($q_{\mathrm{sat}\mathrm{g}}<q_{\mathrm{a}}$), $h_{\mathrm{ug}}=\mathrm{1}$ (see , , for details). The effective field capacity $w_{\mathrm{fc},\mathrm{1}}^{\ast }$ is computed relative to the liquid water content of the uppermost soil layer (it is adjusted in the presence of soil ice compared to the default field capacity). The analogous form holds for the humidity factor over the frozen part of the surface soil layer, $h_{\mathrm{ugf}}$, with $w_{\mathrm{g},\mathrm{1}}$ and $w_{\mathrm{fc},\mathrm{1}}^{\ast }$ replaced by $w_{\mathrm{gf},\mathrm{1}}$ and $w_{\mathrm{fcf},\mathrm{1}}^{\ast }$ (m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}}$) in Eq. (), respectively . Note that it would be more accurate to use $q_{\mathrm{sati}}$ in place of $q_{\mathrm{sat}}$ for the sublimation of the canopy-intercepted snow and the soil ice in Eqs. ()–(), respectively, but this complicates the linearization and this has been neglected for now. The snow factor is defined as $h_{\mathrm{sn}}=L_{\mathrm{s}}/L$. This factor can be modified so that $E_{\mathrm{n}}$ includes both sublimation and evaporation , but the impact of including a liquid water flux has been found to be negligible; thus, for simplicity only sublimation is accounted for currently.

The leading coefficient for the canopy evapotranspiration is defined as $$h_{\mathrm{sv}}=\left(\mathrm{1},-,p_{\mathrm{nv}}\right)h_{\mathrm{svg}}\left(L_{\mathrm{v}},/,L\right)+p_{\mathrm{nv}}{h}_{\mathrm{svn}}\left(L_{\mathrm{s}},/,L\right),$$ where $p_{\mathrm{nv}}$ is an evaporative efficiency factor that is used to partition the canopy interception storage mass flux between evaporation of liquid water and sublimation (see Eq. ). When part of the vegetation canopy is buried (i.e., $p_{\mathrm{n}\mathit{\alpha }}>\mathrm{0}$), a different roughness is felt by the canopy air space so that a new resistance is computed over the $p_{n\mathit{\alpha }}$-covered part of the canopy as is done for sensible heat flux. This is accounted for by defining $$\begin{array}{cc}{\displaystyle }{h}_{\mathrm{svg}}=& {\displaystyle }{p}_{\mathrm{ng}}\left(\mathrm{1},-,p_{n\mathit{\alpha }}\right)\left({\displaystyle \frac{R_{\mathrm{a}\mathrm{v}-\mathrm{c}}}{R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}}}\right)h_{\mathrm{vn}}\\ & {\displaystyle }& {\displaystyle }+\left(\mathrm{1},-,p_{\mathrm{ng}}\right)\left({\displaystyle \frac{R_{\mathrm{a}\mathrm{v}-\mathrm{c}}}{R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}}}\right)h_{\mathrm{vg}},{\displaystyle }{h}_{\mathrm{svn}}=& {\displaystyle }{p}_{\mathrm{ng}}\left(\mathrm{1},-,p_{n\mathit{\alpha }}\right)\left({\displaystyle \frac{R_{\mathrm{a}\mathrm{v}-\mathrm{c}}}{R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}}}\right)\\ & {\displaystyle }& {\displaystyle }+\left(\mathrm{1},-,p_{\mathrm{ng}}\right)\left({\displaystyle \frac{R_{\mathrm{a}\mathrm{v}-\mathrm{c}}}{R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}}}\right).\end{array}$$ The Halstead coefficients in Eq. () are defined as $$\begin{array}{}& {\displaystyle }{h}_{\mathrm{vg}}=& {\displaystyle }\left({\displaystyle \frac{R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}}{R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}+R_{\mathrm{s}}}}\right)\left(\mathrm{1},-,\mathit{\delta }\right)+\mathit{\delta },& {\displaystyle },h_{\mathrm{vn}}=& {\displaystyle }\left({\displaystyle \frac{R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}}{R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}+R_{\mathrm{sn}}}}\right)\left(\mathrm{1},-,\mathit{\delta }\right)+\mathit{\delta }.\end{array}$$ The stomatal resistance $R_{\mathrm{s}}$ can be computed using either the Jarvis method described by or a more physically based method that includes a representation of photosynthesis . The stomatal resistance for the partially snow-buried portion defined as $$\begin{array}{cc}{\displaystyle }{R}_{\mathrm{sn}}=& {\displaystyle }{R}_{\mathrm{s}}/\left[\mathrm{1},,-,,\text{min},\left(p_{n\mathit{\alpha }},,,\mathrm{\hspace{0.17em}\hspace{0.17em}1},,-,,R_{\mathrm{s}},/,R_{\mathrm{s},\max }\right)\right]\\ & {\displaystyle }& {\displaystyle }\left(R_{\mathrm{sn}},\le ,R_{\mathrm{s},\max }\right)\end{array}$$ so that the effect of coverage by the snowpack is to increase the canopy resistance. Note that when the canopy is not partially or fully buried by ground-based snowpack ($p_{n\mathit{\alpha }}=\mathrm{0}$) and does not contain any intercepted snow ($p_{\mathrm{nv}}=\mathrm{0}$), the leading coefficient for the canopy evapotranspiration simplifies to the Halstead coefficient from the composite version of ISBA $$\begin{array}{cc}{\displaystyle }{h}_{\mathrm{sv}}=& {\displaystyle }\left({\displaystyle \frac{R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}}{R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}+R_{\mathrm{s}}}}\right)\left(\mathrm{1},-,\mathit{\delta }\right)+\mathit{\delta }\\ & {\displaystyle }& {\displaystyle }\left(p_{n\mathit{\alpha }},=,\mathrm{0},,,,\text{and},,,,p_{\mathrm{nv}},=,\mathrm{0}\right).\end{array}$$ The fraction of the vegetation covered by water is $\mathit{\delta }$ and is described in Sect. .

The evapotranspiration from the vegetation canopy, $E_{\mathrm{v}}$, is comprised of three components: $$E_{\mathrm{v}}=E_{\mathrm{tr}}+E_{\mathrm{r}}+E_{\mathrm{rn}},$$ where the transpiration, evaporation from the canopy liquid water interception store, and sublimation from the canopy snow interception store are represented by $E_{\mathrm{tr}}$, $E_{\mathrm{r}}$, and $E_{\mathrm{rn}}$, respectively. The expressions for these fluxes are given in Appendix .

Radiative fluxes

The $R_{\mathrm{n}}$ terms in Eqs. ()–() represent the surface net radiation terms (longwave and shortwave components) $$R_{\mathrm{n}x}={\mathrm{SW}}_{\mathrm{net},x}+{\mathrm{LW}}_{\mathrm{net},x},$$ where $x$ $=$ n, g, or v. The total net radiation of the surface is $$\begin{array}{cc}{\displaystyle }{R}_{\mathrm{n}}=& {\displaystyle }{R}_{\mathrm{n}\mathrm{n}}+R_{\mathrm{n}\mathrm{g}}+R_{\mathrm{n}\mathrm{v}}\\ & {\displaystyle }& {\displaystyle }=\mathrm{SW}\downarrow -\mathrm{SW}\uparrow +\mathrm{LW}\downarrow -\mathrm{LW}\uparrow ,\end{array}$$ where the total down-welling solar (shortwave) and atmospheric (longwave) radiative fluxes (W m${}^{-\mathrm{2}}$) at the top of the canopy or snow surface (in the case snow is burying the vegetation) are represented by $\mathrm{SW}\downarrow$ and $\mathrm{LW}\downarrow$, respectively. The total upwelling (towards the atmosphere) shortwave and longwave radiative fluxes $\mathrm{SW}\uparrow$ and $\mathrm{LW}\uparrow$, respectively, are simply defined as the downward components less the total surface net radiative fluxes (summed over the three surfaces). The effective total surface albedo and surface radiative temperature (and emissivity) can then be diagnosed (see the Sect. ) for coupling with the host atmospheric model. The ${\mathit{\tau }}_{\mathrm{n}}$ is defined as the solar radiation transmission at the base of a snowpack layer, so for a sufficiently thin snowpack, solar energy penetrating the snow to the underlying ground surface is expressed as ${\mathit{\tau }}_{\mathrm{n},N_{\mathrm{n}}}{\mathrm{SW}}_{\mathrm{net}\mathrm{n}}$, where $N_{\mathrm{n}}$ represents the number of modeled snowpack layers (for a deep snowpack, this term becomes negligible).

Shortwave radiative fluxes

The total land surface shortwave energy budget can be shown to satisfy $$\mathrm{SW}\downarrow ={\mathrm{SW}}_{\mathrm{net}\mathrm{g}}+{\mathrm{SW}}_{\mathrm{net}\mathrm{v}}+{\mathrm{SW}}_{\mathrm{net}\mathrm{n}}+\mathrm{SW}\uparrow ,$$ where ${\mathrm{SW}}_{\mathrm{net}\mathrm{g}}$, ${\mathrm{SW}}_{\mathrm{net}\mathrm{v}}$, and ${\mathrm{SW}}_{\mathrm{net}\mathrm{n}}$ represent the net shortwave terms for the ground, vegetation canopy, and the ground-based snowpack. The effective surface albedo (which may be required by the atmospheric radiation scheme or for comparison with satellite-based data) is diagnosed as $${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{\mathrm{s}}=\mathrm{SW}\uparrow /\mathrm{SW}\downarrow .$$

The multi-level transmission computations for direct and diffuse radiation are from . The distinction between the visible (VIS) and near-infrared (NIR) radiation components is important in terms of interactions with the vegetation canopy. Here, we take into account two spectral bands for the soil and the vegetation, where visible wavelengths range from approximately 0.3 to 0.7 $\times {10}^{-\mathrm{6}}$ m, and NIR wavelengths range from approximately 0.7 to 1.4 $\times {10}^{-\mathrm{6}}$ m. The spectral values for the soil and the vegetation are provided by ECOCLIMAP as a function of vegetation type and climate.

The effective all-wavelength ground (below-canopy) albedo is defined as $${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{\mathrm{gn}}=p_{\mathrm{ng}}{\mathit{\alpha }}_{\mathrm{n}}+\left(\mathrm{1},-,p_{\mathrm{ng}}\right){\mathit{\alpha }}_{\mathrm{g}},$$ where ${\mathit{\alpha }}_{\mathrm{g}}$ represents the ground albedo.

The ground-based snow albedo, ${\mathit{\alpha }}_{\mathrm{n}}$, is prognostic and depends on the snow grain size. It currently includes up to three spectral bands ; however, when coupled to MEB, only the two aforementioned spectral bands are currently considered for consistency with the vegetation and soil.

The effective canopy albedo, ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{\mathrm{v}}$, represents the combined canopy vegetation, ${\mathit{\alpha }}_{\mathrm{v}}$, and intercepted snow albedos. Currently, however, we assume that ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{\mathrm{v}}={\mathit{\alpha }}_{\mathrm{v}}$, which is based on recommendations by . They showed that multiple reflections and scattering of light from patches of intercepted snow together with a high probability of reflected light reaching the underside of an overlying branch implied that trees actually act like light traps. Thus, they concluded that intercepted snow had no significant influence on the shortwave albedo or the net radiative exchange of boreal conifer canopies.

In addition to baseline albedo values required by the radiative transfer model for each spectral band, the model requires the direct and diffusive downwelling solar components. The diffuse fraction can be provided by observations (offline mode) or a host atmospheric model. For the case when no diffuse information is provided to the surface model, the diffuse fraction is computed using the method proposed by .

Longwave radiative fluxes

The longwave radiation scheme is based on a representation of the vegetation canopy as a plane-parallel surface. The model considers one reflection with three reflecting surfaces (ground, ground-based snowpack, and the vegetation canopy; a schematic is shown in Appendix ). The total land surface longwave energy budget can be shown to satisfy $$\mathrm{LW}\downarrow ={\mathrm{LW}}_{\mathrm{net}\mathrm{g}}+{\mathrm{LW}}_{\mathrm{net}\mathrm{v}}+{\mathrm{LW}}_{\mathrm{net}\mathrm{n}}+\mathrm{LW}\uparrow ,$$ where ${\mathrm{LW}}_{\mathrm{net}\mathrm{g}}$, ${\mathrm{LW}}_{\mathrm{net}\mathrm{v}}$, and ${\mathrm{LW}}_{\mathrm{net}\mathrm{n}}$ represent the net longwave terms for the ground, vegetation canopy, and the ground-based snowpack. The effective surface radiative temperature (which may be required by the atmospheric radiation scheme or for comparison with satellite-based data) is diagnosed as $$T_{\mathrm{rad}}={\left[{\displaystyle \frac{\mathrm{LW}\uparrow -\mathrm{LW}\downarrow \left(\mathrm{1},-,{\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{s}}\right)}{{\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{s}}\mathit{\sigma }}}\right]}^{\mathrm{1}/\mathrm{4}},$$ where $\mathit{\sigma }$ is the Stefan–Boltzmann constant, and ${\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{s}}$ represents the effective surface emissivity. In Eq. (), there are two that are known ($\mathrm{LW}$ fluxes) and two that are unknown ($T_{\mathrm{rad}}$ and ${\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{s}}$). Here we opt to pre-define ${\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{s}}$ in a manner that is consistent with the various surface contributions as $${\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{s}}=p_{\mathrm{ng}}{\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{sn}}+\left(\mathrm{1},-,p_{\mathrm{ng}}\right){\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{sg}}.$$ The canopy-absorption-weighted effective snow and ground emissivities are defined, respectively, as $$\begin{array}{}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{sn}}=& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{n}\mathrm{LW}}{\mathit{ϵ}}_{\mathrm{v}}+\left(\mathrm{1},-,{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{n}\mathrm{LW}}\right){\mathit{ϵ}}_{\mathrm{n}},& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{sg}}=& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{g}\mathrm{LW}}{\mathit{ϵ}}_{\mathrm{v}}+\left(\mathrm{1},-,{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{g}\mathrm{LW}}\right){\mathit{ϵ}}_{\mathrm{g}},\end{array}$$ where ${\mathit{ϵ}}_{\mathrm{v}}$, ${\mathit{ϵ}}_{\mathrm{g}}$, and ${\mathit{ϵ}}_{\mathrm{n}}$ represent the emissivities of the vegetation, snow-free ground, and the ground-based snowpack, respectively. The ground and vegetation emissivities are given by ECOCLIMAP are vary primarily as a function of vegetation class for spatially distributed simulations, or they can be prescribed for local scale studies. The snow emissivity is currently defined as ${\mathit{ϵ}}_{\mathrm{n}}=0.99$. The effect of longwave absorption through the non-snow-buried part of the vegetation canopy is included as $$\begin{array}{cc}{\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{n}\mathrm{LW}}=& {\displaystyle }\left[\mathrm{1},,-,,p_{\mathrm{ng}},,-,,p_{\mathrm{n}\mathit{\alpha }},\left(\mathrm{1},-,p_{\mathrm{ng}}\right)\right]{\mathit{\sigma }}_{\mathrm{LW}}+\left[p_{\mathrm{ng}}\right.\\ & {\displaystyle }& {\displaystyle }\left..,+,,p_{\mathrm{n}\mathit{\alpha }},\left(\mathrm{1},-,p_{\mathrm{ng}}\right)\right]{\mathit{\sigma }}_{\mathrm{f}\mathrm{LW}},& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{g}\mathrm{LW}}=& {\displaystyle }\left[\mathrm{1},,-,,p_{\mathrm{ng}},\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right)\right]{\mathit{\sigma }}_{\mathrm{LW}}+p_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right){\mathit{\sigma }}_{\mathrm{f}\mathrm{LW}},\end{array}$$ where the canopy absorption is defined as $${\mathit{\sigma }}_{\mathrm{LW}}=\mathrm{1}-\exp \left(-,{\mathit{\tau }}_{\mathrm{LW}},,\mathrm{LAI}\right)=\mathrm{1}-{\mathit{\chi }}_{\mathrm{v}}$$ and ${\mathit{\tau }}_{\mathrm{LW}}$ represents a longwave radiation transmission factor that can be species (or land classification) dependent, ${\mathit{\chi }}_{\mathrm{v}}$ is defined as a vegetation view factor, and LAI represents the leaf area index (m${}^{\mathrm{2}}$ m${}^{-\mathrm{2}}$). The absorption over the understory snow-covered fraction of the grid box is modeled quite simply from Eq. () $$\begin{array}{cc}{\displaystyle }{\mathit{\sigma }}_{\mathrm{f}\mathrm{LW}}& {\displaystyle }=\mathrm{\hspace{0.17em}1}-\exp \left[-,{\mathit{\tau }}_{\mathrm{LW}},,\mathrm{LAI},\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right)\right]\\ & {\displaystyle }& {\displaystyle }=\mathrm{\hspace{0.17em}1}-\exp \left[-,{\mathit{\tau }}_{\mathrm{LW}},,{\mathrm{LAI}}_{\mathrm{n}}\right]\end{array}$$ so that transmission is unity (no absorption or reflection by the canopy: ${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{LW}}={\mathit{\sigma }}_{\mathrm{f}\mathrm{LW}}=\mathrm{0}$) when $p_{\mathrm{n}\mathit{\alpha }}=\mathrm{1}$ (i.e., when the canopy has been buried by snow); ${\mathrm{LAI}}_{\mathrm{n}}$ is used to represent the LAI, which has been reduced owing to burial by the snowpack. From Eqs. ()–(), it can be seen that when there is no snowpack (i.e., $p_{\mathrm{ng}}=\mathrm{0}$ and $p_{\mathrm{n}\mathit{\alpha }}=\mathrm{0}$), then the effective surface emissivity is simply an absorption-weighted soil–vegetation value defined as ${\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{s}}={\mathit{\sigma }}_{\mathrm{LW}}{\mathit{ϵ}}_{\mathrm{v}}+\left(\mathrm{1},-,{\mathit{\sigma }}_{\mathrm{LW}}\right){\mathit{ϵ}}_{\mathrm{g}}$. See Appendix  for the derivation of the net longwave radiation terms in Eq. ().

Heat conduction fluxes

The sub-surface snow and ground heat conduction fluxes are modeled using Fourier's law ($G=\mathit{\lambda }\partial T/\partial z$). The heat conduction fluxes in Eqs. ()–() are written in discrete form as $$\begin{array}{}& {\displaystyle }{G}_{\mathrm{g},\mathrm{1}}& {\displaystyle }={\displaystyle \frac{\mathrm{2}\left(T_{\mathrm{g},\mathrm{1}},-,T_{\mathrm{g},\mathrm{2}}\right)}{\left(\mathrm{\Delta },z_{\mathrm{g},\mathrm{1}},/,{\mathit{\lambda }}_{\mathrm{g},\mathrm{1}}\right)+\left(\mathrm{\Delta },z_{\mathrm{g},\mathrm{2}},/,{\mathit{\lambda }}_{\mathrm{g},\mathrm{2}}\right)}}={\mathrm{\Lambda }}_{\mathrm{g},\mathrm{1}}\left(T_{\mathrm{g},\mathrm{1}},-,T_{\mathrm{g},\mathrm{2}}\right),& {\displaystyle }{G}_{\mathrm{n},\mathrm{1}}& {\displaystyle }={\displaystyle \frac{\mathrm{2}\left(T_{\mathrm{n},\mathrm{1}},-,T_{\mathrm{n},\mathrm{2}}\right)}{\left(D_{\mathrm{n},\mathrm{1}},/,{\mathit{\lambda }}_{\mathrm{n},\mathrm{1}}\right)+\left(D_{\mathrm{n},\mathrm{2}},/,{\mathit{\lambda }}_{\mathrm{n},\mathrm{2}}\right)}}={\mathrm{\Lambda }}_{\mathrm{n},\mathrm{1}}\left(T_{\mathrm{n},\mathrm{1}},-,T_{\mathrm{n},\mathrm{2}}\right),& {\displaystyle }{G}_{\mathrm{gn}}& {\displaystyle }={\displaystyle \frac{\mathrm{2}\left(T_{\mathrm{n},N_{\mathrm{n}}},-,T_{\mathrm{g},\mathrm{1}}\right)}{\left(D_{\mathrm{n},N_{\mathrm{n}}},/,{\mathit{\lambda }}_{\mathrm{n},N_{\mathrm{n}}}\right)+\left(\mathrm{\Delta },z_{\mathrm{g},\mathrm{1}},/,{\mathit{\lambda }}_{\mathrm{g},\mathrm{1}}\right)}}={\mathrm{\Lambda }}_{\mathrm{g},\mathrm{n}}\left(T_{\mathrm{n},N_{\mathrm{n}}},-,T_{\mathrm{g},\mathrm{1}}\right),\end{array}$$ where $G_{\mathrm{gn}}$ represents the snow–ground inter-facial heat flux which defines the snow scheme lower boundary condition. All of the internal heat conduction fluxes ($k=\mathrm{2},N-\mathrm{1}$) use the same form as in Eq. () for the snow and Eq. () for the soil . The heat capacities and thermal conductivities ${\mathit{\lambda }}_{\mathrm{g}}$ for the ground depend on the soil texture, organic content , and potentially on the thermal properties of the forest litter in the uppermost layer ; all of the aforementioned properties depend on the water content. The snow thermal property parameterization is described in .

Aerodynamic resistances

The resistances between the surface and the overlying atmosphere, $R_{\mathrm{a}\mathrm{n}-\mathrm{a}}$ and $R_{\mathrm{a}\mathrm{c}-\mathrm{a}}$, are based on modified by to account for different roughness length values for heat and momentum as in ISBA: the full expressions are given in .

Aerodynamic resistance between the bulk vegetation layer and the canopy air

The aerodynamic resistance between the vegetation canopy and the surrounding airspace can be defined as $$R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}={\left(g_{\mathrm{av}},,+,,g_{\mathrm{av}}^{\ast }\right)}^{-\mathrm{1}}.$$ The parameterization of the bulk canopy aerodynamic conductance $g_{\mathrm{av}}$ between the canopy and the canopy air is based on . It is defined as $$g_{\mathrm{av}}={\displaystyle \frac{\mathrm{2}\mathrm{LAI}{a}_{\mathrm{av}}}{{\mathit{\varphi }}_{\mathrm{v}}^{\prime }}}{\left({\displaystyle \frac{u_{\mathrm{hv}}}{\mathrm{lw}}}\right)}^{\mathrm{1}/\mathrm{2}}[\mathrm{1}-\exp (-{\mathit{\varphi }}_{\mathrm{v}}^{\prime }/\mathrm{2}\left)\right],$$ where $u_{\mathrm{hv}}$ represents the wind speed at the top of the canopy (m s${}^{-\mathrm{1}}$),and lw represents the leaf width (m: see Table ). The remaining parameters and their values are defined in Table . The conductance accounting for the free convection correction from is expressed as $$g_{\mathrm{av}}^{\ast }=\left[{\displaystyle \frac{\mathrm{LAI}}{890}},{\left({\displaystyle \frac{T_{\mathrm{v}}-T_{\mathrm{c}}}{\mathrm{lw}}}\right)}^{\mathrm{1}/\mathrm{4}}\right]\left(T_{\mathrm{v}},\ge ,T_{\mathrm{c}}\right).$$ Note that this correction is only used for unstable conditions. The effect of snow burying the vegetation impacts the aerodynamic resistance of the canopy is simply modeled by modifying the LAI using $${\mathrm{LAI}}_{\mathrm{n}}=\mathrm{LAI}\left(\mathrm{1},,-,,p_{\mathrm{n}\mathit{\alpha }}\right).$$ The ${\mathrm{LAI}}_{\mathrm{n}}$ is then used in Eq. () to compute $R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}$, and this resistance is limited to 5000 s m${}^{-\mathrm{1}}$ as ${\mathrm{LAI}}_{\mathrm{n}}\to \mathrm{0}$.

Surface vegetation canopy turbulence parameters that are constant.

Symbol	Definition	Unit	Value	Reference	Comment
$a_{\mathrm{av}}$	Canopy conductance scale factor	m s${}^{-\mathrm{1}/\mathrm{2}}$	0.01		Eq. (26)
${\mathit{\varphi }}_{\mathrm{v}}^{\prime }$	Attenuation coeff. for wind	–	3		p. 386
lw	Leaf width	m	0.02
${\mathit{\varphi }}_{\mathrm{v}}$	Attenuation coeff. for mom.	–	2		p. 386
$z_{\mathrm{0}\mathrm{g}}$	Roughness of soil surface	m	0.007
${\mathit{\chi }}_{\mathrm{L}}$	Ross–Goudriaan leaf angle dist.	–	0.12		p. 26
$u_{\mathrm{l}}$	Typical local wind speed	m s${}^{-\mathrm{1}}$	1		Eq. (B7)
$\mathit{\upsilon }$	Kinematic viscos. of air	m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$	$0.15\times {10}^{-\mathrm{4}}$

Aerodynamic resistance between the ground and the canopy air

The resistance between the ground and the canopy air space is defined as $$R_{\mathrm{ag}-\mathrm{c}}=R_{\mathrm{ag}\mathrm{n}}/{\mathit{\psi }}_H,$$ where $R_{\mathrm{ag}\mathrm{n}}$ is the default resistance value for neutral conditions. The stability correction term ${\mathit{\psi }}_H$ depends on the canopy structural parameters, wind speed, and temperature gradient between the surface and the canopy air. The aerodynamic resistance is also based on . It is assumed that the eddy diffusivity $K$ (m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$) in the vegetation layer follows an exponential profile: $$K\left(z\right)=K\left(z_{\mathrm{hv}}\right)\exp \left[{\mathit{\varphi }}_{\mathrm{v}},\left(\mathrm{1},-,{\displaystyle \frac{z}{z_{\mathrm{hv}}}}\right)\right],$$ where $z_{\mathrm{hv}}$ represents the canopy height. Integrating the reciprocal of the diffusivity defined in Eq. () from $z_{\mathrm{0}\mathrm{g}}$ to $d+z_{\mathrm{0}\mathrm{v}}$ yields $$\begin{array}{cc}{\displaystyle }{R}_{\mathrm{ag}\mathrm{n}}=& {\displaystyle \frac{z_{\mathrm{hv}}}{{\mathit{\varphi }}_{\mathrm{v}}K\left(z_{\mathrm{hv}}\right)}}\mathit{\{}\exp \left[{\mathit{\varphi }}_{\mathrm{v}},\left(\mathrm{1},-,{\displaystyle \frac{z_{\mathrm{0}\mathrm{g}}}{z_{\mathrm{hv}}}}\right)\right]\\ & {\displaystyle }& {\displaystyle }-\exp \left[{\mathit{\varphi }}_{\mathrm{v}},\left(\mathrm{1},-,{\displaystyle \frac{d+z_{\mathrm{0}\mathrm{v}}}{z_{\mathrm{hv}}}}\right)\right]\mathit{\}}.\end{array}$$ The diffusivity at the canopy top is defined as $$K\left(z_{\mathrm{hv}}\right)=k{u}_{\ast \mathrm{hv}}\left(z_{\mathrm{hv}},-,d\right).$$ The von Karman constant $k$ has a value of 0.4. The displacement height is defined as $$d=1.1z_{\mathrm{hv}}\ln \left[\mathrm{1},+,{\left(c_{\mathrm{d}},,\mathrm{LAI}\right)}^{\mathrm{1}/\mathrm{4}}\right],$$ where the leaf drag coefficient $c_{\mathrm{d}}$,is defined from $$c_{\mathrm{d}}=1.328\left[{\displaystyle \frac{\mathrm{2}}{{R_e}^{\mathrm{1}/\mathrm{2}}}}\right]+0.45{\left[{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}},(,\mathrm{1},-,{\mathit{\chi }}_{\mathrm{L}},)\right]}^{1.6},$$ where ${\mathit{\chi }}_{\mathrm{L}}$ represents the Ross–Goudriaan leaf angle distribution function, which has been estimated according to (see Table ), and $R_e$ is the Reynolds number defined as $$R_e={\displaystyle \frac{u_{\mathrm{l}}\mathrm{lw}}{\mathit{\upsilon }}}.$$

The friction velocity at the top of the vegetation canopy is defined as $$u_{\ast \mathrm{hv}}={\displaystyle \frac{k{u}_{\mathrm{hv}}}{\ln \left[\left(z_{\mathrm{hv}},-,d\right),/,z_{\mathrm{0}\mathrm{v}}\right]}},$$ where the wind speed at the top of the canopy is $$u_{\mathrm{hv}}=f_{\mathrm{hv}}{V}_{\mathrm{a}}$$ and $V_{\mathrm{a}}$ represents the wind speed at the reference height $z_{\mathrm{a}}$ above the canopy. The canopy height is defined based on vegetation class and climate within ECOCLIMAP as a primary parameter. It can also be defined using an external dataset, such as from a satellite-derived product (as a function of space and time). The vegetation roughness length for momentum is then computed as a secondary parameter as a function of the vegetation canopy height. The factor $f_{\mathrm{hv}}$ ($\le \mathrm{1}$) is a stability-dependent adjustment factor (see Appendix ).

The dimensionless height scaling factor is defined as $${\mathit{\varphi }}_z={\displaystyle \frac{\left(z_{\mathrm{hv}},-,d\right)}{z_{\mathrm{r}}}}\left({\mathit{\varphi }}_z,\le ,\mathrm{1}\right).$$ The reference height is defined as $z_{\mathrm{r}}=z_{\mathrm{a}}-d$ for simulations where the reference height is sufficiently above the top of the vegetation canopy. This is usually the case for local scale studies using observation data. When MEB is coupled to an atmospheric model, however, the lowest model level can be below the canopy height; therefore, for coupled model simulations $z_{\mathrm{r}}=\max \left(z_{\mathrm{a}},,,,z_{\mathrm{hv}},-,d,+,z_{min}\right)$ where $z_{min}=\mathrm{2}$ (m).

Finally, the stability correction factor from Eq. () is defined as $$\begin{array}{cc}& {\displaystyle }{\mathit{\psi }}_H& {\displaystyle }={\left(\mathrm{1},,-,,a_{\mathrm{hv}},,R_i\right)}^{\mathrm{1}/\mathrm{2}}\left(R_i,\le ,\mathrm{0}\right),{\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+b{R}_i{\left(\mathrm{1},+,c,,R_i\right)}^{\mathrm{1}/\mathrm{2}}}}\left[\mathrm{1},,+,,\left({\displaystyle \frac{R_i}{R_{i,\mathrm{crit}}}}\right),\left(f_{z\mathrm{0}},-,\mathrm{1}\right)\right]\\ & {\displaystyle }& {\displaystyle }\left(R_i,>,\mathrm{0},,,\mathrm{and},,,R_i,\le ,R_{i,\mathrm{crit}}\right),& {\displaystyle }& {\displaystyle }={\displaystyle \frac{f_{z\mathrm{0}}}{\mathrm{1}+b{R}_i{\left(\mathrm{1},+,c,,R_i\right)}^{\mathrm{1}/\mathrm{2}}}}\left(R_i,>,R_{i,\mathrm{crit}}\right),\end{array}$$ where the Richardson number is defined as $$R_i={\displaystyle \frac{-g{z}_{\mathrm{hv}}\left(T_{\mathrm{s}},-,T_{\mathrm{c}}\right)}{T_{\mathrm{s}}{u_{\mathrm{hv}}}^{\mathrm{2}}}}.$$ Note that strictly speaking, the temperature factor in the denominator should be defined as $\left(T_{\mathrm{s}},+,T_{\mathrm{c}}\right)/\mathrm{2}$, but this has only a minor impact for our purposes. The critical Richardson number, $R_{i,\mathrm{crit}}$, is set to 0.2. This parameter has been defined assuming that some turbulent exchange is likely always present (even if intermittent), but it is recognized that eventually a more robust approach should be developed for very stable surface layers . The expression for unstable conditions (Eq. ) is from where the structural parameter is defined as $a_{\mathrm{hv}}=\mathrm{9}$.

It is generally accepted that there is a need to improve the parameterization of the exchange coefficient for extremely stable conditions typically encountered over snow . Since the goal here is not to develop a new parameterization, we simply modify the expression for stable conditions by using the standard function from ISBA. The standard ISBA stability correction for stable conditions is given by Eq. (), where $b=15$ and $c=\mathrm{5}$ . The factor that takes into account differing roughness lengths for heat and momentum is defined as $$f_{z\mathrm{0}}={\displaystyle \frac{\ln \left(z_{\mathrm{hv}},/,z_{\mathrm{0}\mathrm{g}}\right)}{\ln \left(z_{\mathrm{hv}},/,z_{\mathrm{0}\mathrm{gh}}\right)}},$$ where $z_{\mathrm{0}\mathrm{gh}}$ is the ground roughness length for scalars. The weighting function (i.e., ratio of $R_i$ to $R_{i,\mathrm{crit}}$) in Eq. () is used in order to avoid a discontinuity at $R_i=\mathrm{0}$ (the roughness length factor effect vanishes at $R_i=\mathrm{0}$) in Eq. (). An example of Eq.  is shown in Fig.  using the $z_{\mathrm{0}\mathrm{g}}$ from Table , and for $z_{\mathrm{0}\mathrm{gh}}/z_{\mathrm{0}\mathrm{g}}$ of 0.1 and 1.0. Finally, the resistance between the ground-based snowpack $R_{\mathrm{a}\mathrm{n}-\mathrm{c}}$ and the canopy air use the same expressions as for the aerodynamic resistance between the ground and the canopy air outlined herein, but with the surface properties of the snowpack (namely the roughness length and snow surface temperature).

Stability correction term is shown using the Sellers formulation for $R_i\le \mathrm{0}$ while the function for stable conditions adapted from ISBA ($R_i>\mathrm{0}$) for two ratios of $z_{\mathrm{0}\mathrm{g}}/z_{\mathrm{0}\mathrm{gh}}$. The ground surface roughness length is defined in Table .

[Figure omitted. See PDF]

Ground resistance

The soil resistance term is defined based on as $$R_{\mathrm{g}}=\exp \left[a_{\mathrm{Rg}},-,b_{\mathrm{Rg}},,\left({\stackrel{\mathrm{‾}}{w}}_{\mathrm{g}},/,{\stackrel{\mathrm{‾}}{w}}_{\mathrm{sat}}\right)\right].$$ The coefficients are $a_{\mathrm{Rg}}=8.206$ and $b_{\mathrm{Rg}}=4.255$, and the vertically averaged volumetric water content and saturated volumetric water content are given by ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{g}}$ and ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{sat}}$, respectively. The averaging is done from one to several upper layers. Indeed, the inclusion of an explicit ground surface energy budget makes it more conceptually straightforward to include a ground resistance compared to the original composite soil–vegetation surface. The ground resistance is often used as a surrogate for an additional resistance arising due to a forest litter layer, therefore the soil resistance is set to zero when the litter-layer option is activated. Finally, the coefficients $a_{\mathrm{Rg}}$ and $b_{\mathrm{Rg}}$ were determined from a case study for a specific location, and could possibly be location dependent. But currently these values are used, in part, since the litter formulation is the default configuration for MEB for forests as it generally gives better surface fluxes .

Water budget

The governing equations for (water) mass for the bulk canopy, and surface snow and ground layers are written as $$\begin{array}{cc}& {\displaystyle }& {\displaystyle \frac{\partial W_{\mathrm{r}}}{\partial t}}=P_{\mathrm{rv}}+\max \left(\mathrm{0},,,,-,E_{\mathrm{tr}}\right)-E_{\mathrm{r}}-D_{\mathrm{rv}}-{\mathrm{\Phi }}_{\mathrm{v}},& {\displaystyle }& {\displaystyle \frac{\partial W_{\mathrm{r}\mathrm{n}}}{\partial t}}=I_{\mathrm{n}}-U_{\mathrm{n}}-E_{\mathrm{rn}}+{\mathrm{\Phi }}_{\mathrm{v}}{\displaystyle }& {\displaystyle }{p}_{\mathrm{ng}}{\displaystyle \frac{\partial W_{\mathrm{n},\mathrm{1}}}{\partial t}}=P_{\mathrm{s}}-I_{\mathrm{n}}+U_{\mathrm{n}}+p_{\mathrm{ng}}\\ & {\displaystyle }& {\displaystyle }\left(P_{\mathrm{r}},,-,,P_{\mathrm{rv}},,+,,D_{\mathrm{rv}},-,F_{\mathrm{nl},\mathrm{1}},-,E_{\mathrm{n}},+,{\mathrm{\Phi }}_{\mathrm{n},\mathrm{1}},+,{\mathit{\xi }}_{\mathrm{nl},\mathrm{1}}\right),{\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\mathrm{w}}\mathrm{\Delta }{z}_{\mathrm{g},\mathrm{1}}{\displaystyle \frac{\partial w_{\mathrm{g},\mathrm{1}}}{\partial t}}=\left(P_{\mathrm{r}},,-,,P_{\mathrm{rv}},,+,,D_{\mathrm{rv}},,-,,E_{\mathrm{g}}\right)\left(\mathrm{1},-,p_{\mathrm{ng}}\right)\\ & {\displaystyle }& {\displaystyle }+p_{\mathrm{ng}}{F}_{\mathrm{nl},N_{\mathrm{n}}}-R_{\mathrm{0}}-F_{\mathrm{g},\mathrm{1}}-{\mathrm{\Phi }}_{\mathrm{g},\mathrm{1}},& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\mathrm{w}}\mathrm{\Delta }{z}_{\mathrm{g},\mathrm{1}}{\displaystyle \frac{\partial w_{\mathrm{gf},\mathrm{1}}}{\partial t}}={\mathrm{\Phi }}_{\mathrm{g},\mathrm{1}}-E_{\mathrm{gf}}\left(\mathrm{1},-,p_{\mathrm{ng}}\right),\end{array}$$ where $W_{\mathrm{r}}$ and $W_{\mathrm{r}\mathrm{n}}$ represent the vegetation canopy water stores (intercepted water) and the intercepted snow and frozen water (all in kg m${}^{-\mathrm{2}}$), respectively. $W_{\mathrm{n},\mathrm{1}}$ represents the snow liquid water equivalent (SWE) for the uppermost snow layer of the multi-layer scheme. The soil liquid water content and water content equivalent of frozen water are defined as $w_{\mathrm{g}}$ and $w_{\mathrm{gf}}$, respectively (m${}^{\mathrm{3}}$ m${}^{-\mathrm{3}}$).

The interception reservoir $W_{\mathrm{r}}$ is modeled as single-layer bucket, with losses represented by evaporation $E_{\mathrm{r}}$, and canopy drip $D_{\mathrm{rv}}$ of liquid water that exceeds a maximum holding capacity (see Sect.  for details). Sources include condensation (negative $E_{\mathrm{r}}$ and $E_{\mathrm{tr}}$) and $P_{\mathrm{rv}}$, which represents the intercepted precipitation. The positive part of $E_{\mathrm{tr}}$ is extracted from the sub-surface soil layers as a function of soil moisture and a prescribed vertical root zone distribution . This equation is the same as that used in ISBA, except for the addition of the phase change term, ${\mathrm{\Phi }}_{\mathrm{v}}$ (kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$). This term has been introduced owing to the introduction of an explicit canopy snow interception reservoir $W_{\mathrm{r}\mathrm{n}}$; the canopy snow and liquid water reservoirs can exchange mass via this term which is modeled as melt less freezing. The remaining rainfall ($P_{\mathrm{r}}-P_{\mathrm{rv}}$) is partitioned between the snow-free and snow-covered ground surface, where $P_{\mathrm{r}}$ represents the total grid cell rainfall rate. The canopy snow interception is more complex, and represents certain baseline processes such as snow interception $I_{\mathrm{n}}$ and unloading $U_{\mathrm{n}}$; see Sect.  for details.

The soil water and snow liquid water vertical fluxes at the base of the surface ground and snow are represented, respectively, by $F_{\mathrm{g},\mathrm{1}}$ using Darcy's Law and by $F_{\mathrm{nl},\mathrm{1}}$ using a tipping-bucket scheme (kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$). The liquid water flux at the base of the snowpack $F_{\mathrm{nl},N_{\mathrm{n}}}$ is directed downward into the soil and consists in the liquid water in excess of the lowest model liquid water-holding capacity. A description of the snow and soil schemes are given in and , respectively. $R_{\mathrm{0}}$ is the surface runoff. It accounts for sub-grid heterogeneity of precipitation, soil moisture, and for when potential infiltration exceeds a maximum rate . The soil liquid water equivalent ice content can have some losses owing to sublimation in the uppermost soil layer $E_{\mathrm{gf}}$ but it mainly evolves owing to phase changes from soil water freeze–thaw ${\mathrm{\Phi }}_{\mathrm{g}}$. The remaining symbols in Eqs. ()–() are defined and described in Sect.  and .

Precipitation interception

Canopy snow interception

The intercepted snow mass budget is described by Eq. (), while the energy budget is included as a part of the bulk canopy prognostic equation (Eq. ). The positive mass contributions acting to increase intercepted snow on canopy are snowfall interception $I_{\mathrm{n}}$, water on canopy that freezes ${\mathrm{\Phi }}_{\mathrm{v}}<\mathrm{0}$, and sublimation of water vapor to ice $E_{\mathrm{rn}}<\mathrm{0}$. Unloading $U_{\mathrm{n}}$, sublimation $E_{\mathrm{rn}}>\mathrm{0}$, and snowmelt ${\mathrm{\Phi }}_{\mathrm{v}}>\mathrm{0}$, are the sinks. All of the terms are in kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$. It is assumed that intercepted rain and snow can co-exist on the canopy. The intercepted snow is assumed to have the same temperature as the canopy $T_{\mathrm{v}}$; thus, there is no advective heat exchange with the atmosphere that simplifies the equations. For simplicity, when intercepted water on the canopy freezes, it is assumed to become part of the intercepted snow.

The parameterization of interception efficiency is based upon . It determines how much snow is intercepted during the time step and is defined as $$I_{\mathrm{n},\mathrm{v},\mathrm{0}}=\left(W_{\mathrm{r}\mathrm{n}}^{\ast },-,W_{\mathrm{r}\mathrm{n}}\right)\left[\mathrm{1},-,\exp ,\left(-,k_{\mathrm{n},\mathrm{v}},,P_{\mathrm{s}},,\mathrm{\Delta },t\right)\right],$$ where ${W_{\mathrm{r}\mathrm{n}}}^{\ast }$ is the maximum snow load allowed, $P_{\mathrm{s}}$ the frozen precipitation rate, and $k_{\mathrm{n},\mathrm{v}}$ a proportionality factor. $k_{\mathrm{n},\mathrm{v}}$ is a function of ${W_{\mathrm{r}\mathrm{n}}}^{\ast }$ and the maximum plan area of the snow–leaf contact area per unit area of ground $C_{\mathrm{n},\mathrm{vp}}$: $$k_{\mathrm{n},\mathrm{v}}={\displaystyle \frac{C_{\mathrm{n},\mathrm{vp}}}{{W_{\mathrm{r}\mathrm{n}}}^{\ast }}}.$$

For a closed canopy, $C_{\mathrm{n},\mathrm{vp}}$ would be equal to one, but for a partly open canopy it is described by the relationship $$C_{\mathrm{n},\mathrm{vp}}={\displaystyle \frac{C_{\mathrm{n},\mathrm{vc}}}{\mathrm{1}-C_{\mathrm{c}}{u}_{\mathrm{hv}}{z}_{\mathrm{hv}}/\left(w_{\mathrm{n}},,J_{\mathrm{n}}\right)}},$$ where $C_{\mathrm{n},\mathrm{vc}}$ is the canopy coverage per unit area of ground which can be expressed as $\mathrm{1}-{\mathit{\chi }}_{\mathrm{v}}$ where ${\mathit{\chi }}_{\mathrm{v}}$ is the sky-view factor (see Eq. ), and $u_{\mathrm{hv}}$ represents the mean horizontal wind speed at the canopy top (Eq. ), which corresponds to the height $z_{\mathrm{hv}}$ (m). The characteristic vertical snow-flake velocity, $w_{\mathrm{n}}$, is set to 0.8 m s${}^{-\mathrm{1}}$ . $J_{\mathrm{n}}$ is set to 10${}^{\mathrm{3}}$ m, which is assumed to represent the typical size of the mean forested down wind distance.

For calm conditions and completely vertically falling snowflakes, $C_{\mathrm{n},\mathrm{vp}}=C_{\mathrm{c}}$. For any existing wind, snow could be intercepted by the surrounding trees so that high wind speed increases interception efficiency. Generally for open boreal conifer canopies, $C_{\mathrm{n},\mathrm{vc}}<C_{\mathrm{n},\mathrm{vp}}<\mathrm{1}$. Under normal wind speed conditions (i.e., wind speeds larger than 1 m s${}^{-\mathrm{1}}$), $C_{\mathrm{n},\mathrm{vc}}$ (and $C_{\mathrm{n},\mathrm{vp}}$) values are usually close to unity.

The maximum allowed canopy snow load, ${W_{\mathrm{r}\mathrm{n}}}^{\ast }$, is a function of the maximum snow load per unit branch area, $S_{\mathrm{n},\mathrm{v}}$ (kg m${}^{-\mathrm{2}}$), and the leaf area index: $${W_{\mathrm{r}\mathrm{n}}}^{\ast }=S_{\mathrm{n},\mathrm{v}}\mathrm{LAI}$$ where $S_{\mathrm{n},\mathrm{v}}$ is defined as $$S_{\mathrm{n},\mathrm{v}}=\stackrel{\mathrm{‾}}{S_{\mathrm{n},\mathrm{v}}}\left(0.27,+,{\displaystyle \frac{46}{{\mathit{\rho }}_{\mathrm{n},\mathrm{v}}}}\right).$$ Based on measurements, estimated average values for $\stackrel{\mathrm{‾}}{S_{\mathrm{n},\mathrm{v}}}$ of 6.6 for pine and 5.9 kg m${}^{-\mathrm{2}}$ for spruce trees. Because the average value for this parameter only varies by about 10 % across these two fairly common tree species, and ECOCLIMAP does not currently make a clear distinction between these two forest classes, we currently use 6.3 as the default value for all forest classes. ${\mathit{\rho }}_{\mathrm{n},\mathrm{v}}$ is the canopy snow density (kg m${}^{-\mathrm{3}}$) defined by the relationship $${\mathit{\rho }}_{\mathrm{n},\mathrm{v}}=67.92+51.25\exp \left[\left(T_{\mathrm{c}},-,T_{\mathrm{f}}\right),/,2.59\right]\left(T_{\mathrm{c}},\le ,T_{\mathrm{c}\max }\right),$$ where $T_{\mathrm{c}}$ is the canopy air temperature and $T_{\mathrm{c}\max }$ is the temperature corresponding to the maximum snow density. Assuming a maximum snow density of 750 kg m${}^{-\mathrm{3}}$ and solving Eq. () for canopy temperature yields $T_{\mathrm{c}\max }=279.854$ K. This gives values of $S_{\mathrm{n},\mathrm{v}}$ in the range 4–6 kg m${}^{-\mathrm{2}}$.

The water vapor flux between the intercepted canopy snow and the canopy air $E_{\mathrm{rn}}$ (Eq. ), includes the evaporative efficiency $p_{\mathrm{nv}}$. This effect was first described by . In the ISBA-MEB parameterization, the formulation is slightly modified so that it approaches zero when there is no intercepted snow load: $$p_{\mathrm{nv}}={\displaystyle \frac{0.89{S_{\mathrm{nv}}}^{0.3}}{\mathrm{1}+\exp \left[-,4.7,(,S_{\mathrm{nv}},-,0.45,)\right]}},$$ where $S_{\mathrm{nv}}$ is the ratio of snow-covered area on the canopy to the total canopy area $$S_{\mathrm{nv}}={\displaystyle \frac{W_{\mathrm{r}\mathrm{n}}}{{W_{\mathrm{r}\mathrm{n}}}^{\ast }}}\left(\mathrm{0},\le ,S_{\mathrm{nv}},\le ,\mathrm{1}\right).$$ A numerical test is performed to determine if the canopy snow becomes less than zero within one time step due to sublimation. If this is true, then the required mass is removed from the underlying snowpack so that the intercepted snow becomes exactly zero during the time step to ensure a high degree of mass conservation. Note that this adjustment is generally negligible.

The intercepted snow unloading, due to processes such as wind and branch bending, has to be estimated. suggested an experimentally verified exponential decay in load over time $t$, which is used in the parameterization: $$U_{\mathrm{n},\mathrm{v}}=I_{\mathrm{n},\mathrm{v},\mathrm{0}}\exp (-U_{\mathrm{nL}}t)=I_{\mathrm{n},\mathrm{v},\mathrm{0}}{c}_{\mathrm{nL}},$$ where $U_{\mathrm{nL}}$ is an unloading rate coefficient (s${}^{-\mathrm{1}}$) and $c_{\mathrm{nL}}$ the dimensionless unloading coefficient. found that $c_{\mathrm{nL}}=0.678$ was a good approximation that, with a time step of 15 min, gives $U_{\mathrm{nL}}=-4.498\times {10}^{-\mathrm{6}}$ s${}^{-\mathrm{1}}$. A tuned value for the RCA-LSM from the Snow Model Intercomparison Project phase 2 (SnowMIP2) experiments is $U_{\mathrm{nL}}=-3.4254\times {10}^{-\mathrm{6}}$ s${}^{-\mathrm{1}}$, which has been adopted for MEB for now. All unloaded snow is assumed to fall to the ground where it is added to the snow storage on forest ground. Further, corrections to compensate for changes in the original LSM due to this new parameterization have been made for heat capacity, latent heat of vaporization, evapotranspiration, snow storages, and fluxes of latent heat.

Finally, canopy snow will partly melt if the temperature rises above the melting point and become intercepted water, where the intercepted (liquid and frozen) water phase change is simply proportional to the temperature: $${\mathrm{\Phi }}_{\mathrm{v}}={\displaystyle \frac{C_{\mathrm{i}}{W}_{\mathrm{r}\mathrm{n}}}{L_{\mathrm{f}}{\mathit{\tau }}_{\mathrm{\Phi }}}}\left(T_{\mathrm{f}},-,T_{\mathrm{v}}\right)={\displaystyle \frac{C_{\mathrm{i}}{S}_{\mathrm{nv}}{W}_{\mathrm{r}\mathrm{n}}^{\ast }}{L_{\mathrm{f}}{\mathit{\tau }}_{\mathrm{\Phi }}}}\left(T_{\mathrm{f}},-,T_{\mathrm{v}}\right),$$ where ${\mathrm{\Phi }}_{\mathrm{v}}<\mathrm{0}$ signifies melting. $T_{\mathrm{f}}$ represents the melting point temperature (273.15 K) and the characteristic phase change timescale is ${\mathit{\tau }}_{\mathrm{\Phi }}$ (s). If it is assumed that the available heating during the time step for phase change is proportional to canopy biomass via the $\mathrm{LAI}$ then Eq. () can be written (for both melt and refreezing) as $${\mathrm{\Phi }}_{\mathrm{v}}=S_{\mathrm{nv}}{k}_{\mathrm{\Phi }\mathrm{v}}\left(T_{\mathrm{f}},-,T_{\mathrm{v}}\right).$$ Note that if energy is available for melting, the phase change rate is limited by the amount of intercepted snow, and likewise freezing is limited by the amount of intercepted liquid water. The melting of intercepted snow within the canopy can be quite complex, thus currently the simple approach in Eq. () adopted herein. The phase change coefficient was tuned to a value of $k_{\mathrm{\Phi }\mathrm{v}}=5.56\times {10}^{-\mathrm{6}}$ kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}}$ for the SNOWMIP2 experiments with the RCA-LSM. Currently, this value is the default for ISBA-MEB.

Canopy rain interception

The rain intercepted by the vegetation is available for potential evaporation, which means that it has a strong influence on the fluxes of heat and consequently also on the surface temperature. The rate of change of intercepted water on vegetation canopy is described by Eq. (). The rate that water is intercepted by the overstory (which is not buried by the ground-based snow) is defined as $$P_{\mathrm{rv}}=P_{\mathrm{r}}\left(\mathrm{1},-,{\mathit{\chi }}_{\mathrm{v}}\right)\left(\mathrm{1},-,p_{\mathrm{ng}},p_{\mathit{\alpha }\mathrm{n}}\right),$$ where ${\mathit{\chi }}_{\mathrm{v}}$ is a view factor indicating how much of the precipitation that should fall directly to the ground (see Eq. ). The overstory canopy drip rate $D_{\mathrm{rv}}$ is defined simply as the value of water in the reservoir which exceeds the maximum holding capacity: $$D_{\mathrm{rv}}=\max \left(\mathrm{0},,,,W_{\mathrm{rv}},-,W_{\mathrm{rv},\max }\right)/\mathrm{\Delta }t,$$ where the maximum liquid water-holding capacity is defined simply as $$W_{\mathrm{rv},\max }=c_{\mathrm{wrv}}\mathrm{LAI}.$$

Generally speaking, $c_{\mathrm{wrv}}=0.2$ , although it can be modified slightly for certain vegetation cover. Note that Eq. () is first evaluated with $D_{\mathrm{rv}}=\mathrm{0}$, and then the canopy drip is computed as a residual. Thus, the final water amount is corrected by removing the canopy drip or throughfall. This water can then become a liquid water source for the soil and the ground-based snowpack.

The fraction of the vegetation covered with water is defined as $$\begin{array}{cc}{\displaystyle }{\mathit{\delta }}_{\mathrm{v}}& {\displaystyle }=\left(\mathrm{1},-,{\mathit{\omega }}_{\mathrm{rv}}\right){\left({\displaystyle \frac{W_{\mathrm{r}}}{W_{\mathrm{r},\max }}}\right)}^{\mathrm{2}/\mathrm{3}}\\ & {\displaystyle }& {\displaystyle }+{\displaystyle \frac{{\mathit{\omega }}_{\mathrm{rv}}{W}_{\mathrm{r}}}{\left(\mathrm{1},+,a_{\mathrm{rv}},,\mathrm{LAI}\right)W_{\mathrm{r},\max }-a_{\mathrm{rv}}{W}_{\mathrm{r}}}}.\end{array}$$ used the first term on the RHS of Eq. () for relatively low vegetation and the second term for tall vegetation . Currently in ISBA, a weighting function is used which introduces the vegetation height dependence using the roughness length as a proxy from $${\mathit{\omega }}_{\mathrm{rv}}=\mathrm{2}{z}_{\mathrm{0}\mathrm{v}}-\mathrm{\hspace{0.17em}1}\left(\mathrm{0},\le ,{\mathit{\omega }}_{\mathrm{rv}},\le ,\mathrm{1}\right),$$ where the current value for the dimensionless coefficient is $a_{\mathrm{rv}}=\mathrm{2}$.

Halstead coefficient

In the case of wet vegetation, the total plant evapotranspiration is partitioned between the evaporation of intercepted water, and transpiration via stomata by the Halstead coefficient. In MEB, two such coefficients are used for the non-snow-buried and buried parts of the vegetation canopy, $h_{\mathrm{vg}}$ and $h_{\mathrm{vn}}$ (Eqs.  and , respectively). In MEB, the general form of the Halstead coefficient, as defined in , is modified by introducing the factor $k_{\mathrm{v}}$ to take into account the fact that saturated vegetation can transpire, i.e., when ${\mathit{\delta }}_{\mathrm{v}}=\mathrm{1}$ . Thus, for MEB we define $\mathit{\delta }=k_{\mathrm{v}}{\mathit{\delta }}_{\mathrm{v}}$. The intercepted water forms full spheres just touching the vegetation surface when $k_{\mathrm{v}}=\mathrm{0}$, which allows for full transpiration from the whole leaf surface. In contrast, $k_{\mathrm{v}}=\mathrm{1}$ would represent a situation where a water film covers the vegetation completely and no transpiration is allowed. To adhere to the interception model as described above, where the intercepted water exists as droplets, we set the value of $k_{\mathrm{v}}$ to $0.25$. Note that in the case of condensation, i.e., $E<\mathrm{0}$, $h_{\mathrm{v}}=\mathrm{1}$.

Without a limitation of $h_{\mathrm{vg}}$ and $h_{\mathrm{vn}}$, the evaporative demand could exceed the available intercepted water during a time step, especially for the canopy vegetation which experiences a relatively low aerodynamic resistance. To avoid such a situation, a maximum value of the Halstead coefficient is imposed by calculating a maximum value of the ${\mathit{\delta }}_{\mathrm{v}}$. See Appendix  for details.

Conclusions

This paper presents the description of a new multi-energy balance (MEB) scheme for representing tall vegetation in the ISBA land surface model component of the SURFEX land–atmosphere coupling and driving platform. This effort is part of the ongoing effort within the international scientific community to continually improve the representation of land surface processes for hydrological and meteorological research and applications.

MEB consists in a fully implicit numerical coupling between a multi-layer physically based snowpack model, a variable-layer soil scheme, an explicit litter layer, a bulk vegetation scheme, and the atmosphere. It also includes a feature that permits a coupling transition of the snowpack from the canopy air to the free atmosphere as a function of snow depth and canopy height using a fully implicit numerical scheme. MEB has been developed in order to meet the criteria associated with computational efficiency, high coding standards (especially in terms of modularity), conservation (of mass, energy, and momentum), numerical stability for large (time step) scale applications, and state-of-the-art representation of the key land surface processes required for current hydrological and meteorological modeling research and operational applications at Météo-France and within the international community as a part of the HIRLAM consortium. This includes regional scale real-time hindcast hydro-meteorological modeling, coupling within both research and operational non-hydrostatic models, regional climate models, and a global climate model, not to mention being used for ongoing offline land surface reanalysis projects and fundamental research applications.

The simple composite soil–vegetation surface energy budget approach of ISBA has proven its ability to provide solid scientific results and realistic boundary conditions for hydrological and meteorological models since its creation over 2 decades ago. However, owing to the ever increasing demands of the user community, it was decided to improve the representation of the vegetation processes as a priority. The key motivation of the MEB development was to move away from the composite scheme in order to address certain known issues (such as excessive bare-soil evaporation in forested areas, the neglect of canopy snow interception processes), to improve consistency in terms of the representation of the carbon cycle (by modeling explicit vegetation energy and carbon exchanges), to add new key explicit processes (forest litter, the gradual covering of vegetation by ground-based snow cover), and to open the door to potential improvements in land data assimilation (by representing distinct surface temperatures for soil and vegetation). Finally, note that while some LSMs intended for GCMs now use multiple-vegetation layers, a single bulk vegetation layer is currently used in MEB since it has been considered as a reasonable first increase in complexity level from the composite soil–vegetation scheme. However, MEB has been designed such that the addition of more canopy layers could be added if deemed necessary in the future.

This is part one of two companion papers describing the model formulation of ISBA-MEB. Part two describes the model evaluation at the local scale for several contrasting well-instrumented sites in France, and for over 42 sites encompassing a wide range of climate conditions for several different forest classes over multiple annual cycles (Napoly et al., 2016). This two-part series of papers will be followed by a series of papers in upcoming years that will present the evaluation and analysis of ISBA-MEB with a specific focus (coupling with snow processes, regional to global scale hydrology, and finally fully coupled runs in a climate model).

Code availability

The MEB code is a part of the ISBA LSM and is available as open source via the surface modeling platform called SURFEX, which can be downloaded at http://www.cnrm-game-meteo.fr/surfex/. SURFEX is updated at a relatively low frequency (every 3 to 6 months) and the developments presented in this paper are available starting with SURFEX version 8.0. If more frequent updates are needed, or if what is required is not in Open-SURFEX (DrHOOK, FA/LFI formats, GAUSSIAN grid), you are invited to follow the procedure to get a SVN account and to access real-time modifications of the code (see the instructions at the previous link).

Thermodynamic coupling variable

If potential temperature is used as the thermodynamic variable in the coupled model diffusion scheme, then the thermodynamic variable $\mathcal{T}$ (J kg${}^{-\mathrm{1}}$; see Eqs. –) coefficients are defined as $$\begin{array}{}& {\displaystyle }{\mathcal{B}}_x=& {\displaystyle }\mathrm{0}\left(x,=,\mathrm{v},,,\mathrm{g},,,\mathrm{n},,,\mathrm{c},,,\mathrm{a}\right),& {\displaystyle }{\mathcal{A}}_x=& {\displaystyle }{C}_p/{\mathrm{\Pi }}_{\mathrm{s}}\left(x,=,\mathrm{v},,,\mathrm{g},,,\mathrm{n},,,\mathrm{c}\right),& {\displaystyle }{\mathcal{A}}_{\mathrm{a}}=& {\displaystyle }{C}_p/{\mathrm{\Pi }}_{\mathrm{a}},\end{array}$$ where $\mathrm{\Pi }$ is the non-dimensional Exner function and $C_p$ is the heat capacity of dry air (J kg${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}}$). If the atmospheric variable being diffused is dry static energy then $$\begin{array}{}& {\displaystyle }{\mathcal{B}}_x=& {\displaystyle }\mathrm{0}\left(x,=,\mathrm{v},,,\mathrm{g},,,\mathrm{n},,,\mathrm{c}\right),& {\displaystyle }{\mathcal{B}}_{\mathrm{a}}=& {\displaystyle }g{z}_{\mathrm{a}},& {\displaystyle }{\mathcal{A}}_x=& {\displaystyle }{C}_p\left(x,=,v,,,g,,,n,,,c,,,a\right),\end{array}$$ where $z_{\mathrm{a}}$ is the height (m) of the simulated or observed overlying atmospheric temperature, $T_{\mathrm{a}}$ and $g$ is the gravitational constant. The choice of the atmospheric thermodynamic variable is transparent to ISBA-MEB (it is made within the surface–atmosphere coupler). The default (in offline mode and in inline mode with certain atmospheric models) is using Eqs. ()–(). Note that the method can be extended to use the actual air heat capacity (including water vapor) if a linearization of the heat capacity is used.

Latent heat normalization factor

The $L$ is a normalization factor ($L_{\mathrm{v}}\le L\le L_{\mathrm{s}}$), which could be determined in a number of ways. This coefficient ensures conservation of mass between the different surfaces and the atmosphere.

One possible method is to diagnose it by inverting the equation for ${\mathrm{LE}}_{\mathrm{c}}$ (multiplying Eq.  by $L$ thereby eliminating it from the RHS of this equation, and then solving for $L$), but the resulting equation is difficult to apply since the terms can be either positive or negative, and division by a small number is possible. Here, a more smooth (in time) function is proposed, which accounts for each of the surfaces weighted by its respective fraction: $$L={\displaystyle \frac{a_{\mathrm{Ls}}{L}_{\mathrm{s}}+a_{\mathrm{Lv}}{L}_{\mathrm{v}}}{a_{\mathrm{Ls}}+a_{\mathrm{Lv}}}},$$ where $$\begin{array}{cc}& {\displaystyle }{a}_{L_{\mathrm{v}}}& {\displaystyle }=\left[{\mathit{\sigma }}_{\mathrm{f}},,\left(\mathrm{1},-,p_{\mathrm{nv}}\right),,+,,\left(\mathrm{1},-,p_{\mathrm{ng}}\right),\left(\mathrm{1},-,p_{\mathrm{gf}}\right)\right]\left(\mathrm{1},-,p_{\mathrm{ng}},p_{\mathrm{n}\mathit{\alpha }}\right),{\displaystyle }{a}_{L_{\mathrm{s}}}& {\displaystyle }=\left[{\mathit{\sigma }}_{\mathrm{f}},,p_{\mathrm{nv}},,+,,\left(\mathrm{1},-,p_{\mathrm{ng}}\right),,p_{\mathrm{gf}},,+,,p_{\mathrm{ng}}\right]\left(\mathrm{1},-,p_{\mathrm{ng}},p_{\mathrm{n}\mathit{\alpha }}\right)\\ & {\displaystyle }& {\displaystyle }+p_{\mathrm{ng}}{p}_{\mathrm{n}\mathit{\alpha }}.\end{array}$$ In the limit as the snow totally buries the canopy vegetation, $L\to L_{\mathrm{s}}$. In contrast, for snow and surface ice-free conditions, $L=L_{\mathrm{v}}$.

Turbulent flux expressions

The turbulent fluxes of heat and water vapor can be further decomposed into different components, which are required for computing different diagnostics and coupling with the water budgets. They are presented herein.

Sensible heat flux

It is convenient to split $H_{\mathrm{n}}$ into two components since one governs the coupling between the canopy air space and the snow surface, while the other modulates the exchanges with the overlying atmosphere (as the canopy layer becomes buried).

The ground-based snowpack heat flux, $H_{\mathrm{n}}$ (Eq. ), can be split into a part that modulates the heat exchange with the canopy air space, $H_{\mathrm{n}-\mathrm{c}}$ and the other part which controls the exchanges directly with the overlying atmosphere, $H_{\mathrm{n}-\mathrm{a}}$, defined as $$\begin{array}{}& {\displaystyle }{H}_{\mathrm{n}-\mathrm{c}}=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}{\displaystyle \frac{\left({\mathcal{T}}_{\mathrm{n}},-,{\mathcal{T}}_{\mathrm{c}}\right)}{R_{\mathrm{a}\mathrm{n}-\mathrm{c}}}},& {\displaystyle }{H}_{\mathrm{n}-\mathrm{a}}=& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}{\displaystyle \frac{\left({\mathcal{T}}_{\mathrm{n}},-,{\mathcal{T}}_{\mathrm{a}}\right)}{R_{\mathrm{a}\mathrm{n}-\mathrm{a}}}}.\end{array}$$ ${\mathcal{T}}_{\mathrm{c}}$ is diagnosed by imposing conservation of the heat fluxes between the surface and the canopy air (as described in Appendix ). Using the definition in Eq. (), the total sensible heat flux exchange with the atmosphere (Eq. ) can also be written in more compact form as $$H={\mathit{\rho }}_{\mathrm{a}}\left[\left(\mathrm{1},-,p_{\mathrm{ng}},,p_{\mathrm{n}\mathit{\alpha }}\right),,H_{\mathrm{c}},,+,,p_{\mathrm{ng}},,p_{\mathrm{n}\mathit{\alpha }},,H_{\mathrm{n}-\mathrm{a}}\right].$$

Water vapor flux

The various water vapor flux terms must be broken into different components for use within the different water balance equations for the vegetation, soil, and snowpack. Using the definitions in Eqs. ()–(), the components of the canopy evapotranspiration $E_{\mathrm{v}}$ can be expressed as

$$\begin{array}{cc}{\displaystyle }{E}_{\mathrm{tr}}& {\displaystyle }={\mathit{\rho }}_{\mathrm{a}}\left({\displaystyle \frac{L_{\mathrm{v}}}{L}}\right)\left(q_{\mathrm{sat}\mathrm{v}},-,q_{\mathrm{c}}\right)\\ & {\displaystyle }& {\displaystyle }\left[{\displaystyle \frac{p_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right)}{R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}+R_{\mathrm{sn}}}},,+,,{\displaystyle \frac{\mathrm{1}-p_{\mathrm{ng}}}{R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}+R_{\mathrm{s}}}}\right]\left(\mathrm{1},-,p_{\mathrm{nv}}\right)\left(\mathrm{1},-,\mathit{\delta }\right),{\displaystyle }{E}_{\mathrm{r}}& {\displaystyle }={\mathit{\rho }}_{\mathrm{a}}\left({\displaystyle \frac{L_{\mathrm{v}}}{L}}\right)\left(q_{\mathrm{sat}\mathrm{v}},-,q_{\mathrm{c}}\right)\\ & {\displaystyle }& {\displaystyle }\left[{\displaystyle \frac{p_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right)}{R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}}},,+,,{\displaystyle \frac{\mathrm{1}-p_{\mathrm{ng}}}{R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}}}\right]\left(\mathrm{1},-,p_{\mathrm{nv}}\right)\mathit{\delta },{\displaystyle }{E}_{\mathrm{rn}}& {\displaystyle }={\mathit{\rho }}_{\mathrm{a}}\left({\displaystyle \frac{L_{\mathrm{s}}}{L}}\right)\left(q_{\mathrm{sat}\mathrm{v}},-,q_{\mathrm{c}}\right)\\ & {\displaystyle }& {\displaystyle }\left[{\displaystyle \frac{p_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right)}{R_{\mathrm{a}\mathrm{vn}-\mathrm{c}}}},,+,,{\displaystyle \frac{\mathrm{1}-p_{\mathrm{ng}}}{R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}}}\right]p_{\mathrm{nv}}.\end{array}$$ The complex resistances (bracketed terms in Eqs. –) arise owing to the inclusion of the effects of burying the snow canopy by the ground-based snowpack. If the ground-based snowpack is not sufficiently deep to bury any of the canopy ($p_{\mathrm{n}\mathit{\alpha }}=\mathrm{0}$), then the bracketed term in Eq. () simplifies to $\mathrm{1}/\left(R_{\mathrm{a}\mathrm{vg}-\mathrm{c}},+,R_{\mathrm{s}}\right)$ (note that $R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}=R_{\mathrm{a}\mathrm{v}-\mathrm{c}}$ when $p_{\mathrm{n}\mathit{\alpha }}=\mathrm{0}$ from Eq. ), and likewise the bracketed terms in Eqs. ()–() simplify to $\mathrm{1}/R_{\mathrm{a}\mathrm{vg}-\mathrm{c}}$. Finally, the partitioning between the vapor fluxes from intercepted snow and the snow-free canopy reservoir and transpiration is done using $p_{\mathrm{nv}}$, which represents the fraction of the snow interception reservoir that is filled (see Eq. ).

Using the definitions of $q_{\mathrm{g}}$ from Eq. () together with those for the humidity factors, $h_{\mathrm{sg}}$ and $h_{\mathrm{a}}$ (Eqs.  and , respectively) and the soil coefficient, ${\mathit{\delta }}_{\mathrm{g}}$ (Eq. ), the bare-soil evaporation $E_{\mathrm{g}}$ components can be expressed as $$\begin{array}{}& {\displaystyle }{E}_{\mathrm{gl}}& {\displaystyle }={\mathit{\rho }}_{\mathrm{a}}\left({\displaystyle \frac{L_{\mathrm{v}}}{L}}\right)\left(h_{\mathrm{ug}},,q_{\mathrm{sat}\mathrm{g}},-,q_{\mathrm{c}}\right)\left({\displaystyle \frac{{\mathit{\delta }}_{\mathrm{gcor}}}{R_{\mathrm{a}\mathrm{g}}+R_{\mathrm{g}}}}\right)\left(\mathrm{1},-,p_{\mathrm{gf}}\right),& {\displaystyle }{E}_{\mathrm{gf}}& {\displaystyle }={\mathit{\rho }}_{\mathrm{a}}\left({\displaystyle \frac{L_{\mathrm{s}}}{L}}\right)\left(h_{\mathrm{ugf}},,q_{\mathrm{sat}\mathrm{g}},-,q_{\mathrm{c}}\right)\left({\displaystyle \frac{{\mathit{\delta }}_{\mathrm{gfcor}}}{R_{\mathrm{a}\mathrm{g}}+R_{\mathrm{gf}}}}\right)p_{\mathrm{gf}},\end{array}$$ where $E_{\mathrm{g}}=E_{\mathrm{gl}}+E_{\mathrm{gf}}$. The delta function, ${\mathit{\delta }}_{\mathrm{gfcor}}$, is a numerical correction term, which is required owing to the linearization of $q_{\mathrm{sat}\mathrm{g}}$ and is unity unless both $h_{\mathrm{ugf}}{q}_{\mathrm{sat}\mathrm{g}}<q_{\mathrm{c}}$ and $q_{\mathrm{sat}\mathrm{g}}>q_{\mathrm{c}}$, in which case it is set to zero. Note that the ground resistances, $R_{\mathrm{g}}$ and $R_{\mathrm{gf}}$, are set to zero if the forest litter option is active (the default for forests).

The ground-based snowpack sublimation, $E_{\mathrm{n}}$ (Eq. ), can be partitioned into a vapor exchange with the canopy air space, $E_{\mathrm{n}-\mathrm{c}}$ and the overlying atmosphere, $E_{\mathrm{n}-\mathrm{a}}$, as $$\begin{array}{}& {\displaystyle }{E}_{\mathrm{n}-\mathrm{c}}& {\displaystyle }={\mathit{\rho }}_{\mathrm{a}}\left({\displaystyle \frac{L_{\mathrm{s}}}{L}}\right)\left({\displaystyle \frac{q_{\mathrm{sati}\mathrm{n}}-q_{\mathrm{c}}}{R_{\mathrm{a}\mathrm{n}-\mathrm{c}}}}\right),& {\displaystyle }{E}_{\mathrm{n}-\mathrm{a}}& {\displaystyle }={\mathit{\rho }}_{\mathrm{a}}\left({\displaystyle \frac{L_{\mathrm{s}}}{L}}\right)\left({\displaystyle \frac{q_{\mathrm{sati}\mathrm{n}}-q_{\mathrm{a}}}{R_{\mathrm{a}\mathrm{n}-\mathrm{a}}}}\right).\end{array}$$

The corresponding latent heat fluxes can be determined by simply multiplying Eqs. ()–() by $L$. Finally, using the definition in Eq. (), the total vapor exchange with the atmosphere (Eq. ) can also be written in more compact form as $$E={\mathit{\rho }}_{\mathrm{a}}\left[\left(\mathrm{1},-,p_{\mathrm{ng}},,p_{\mathrm{n}\mathit{\alpha }}\right),,E_{\mathrm{c}},,+,,p_{\mathrm{ng}},,p_{\mathrm{n}\mathit{\alpha }},,E_{\mathrm{n}-\mathrm{a}}\right].$$

Canopy-top wind stability factor

The expressions for the stability factor $f_{\mathrm{hv}}$ (Eq. ), which is used to compute the wind at the top of the vegetation canopy $u_{\mathrm{hv}}$, are taken from . They are defined as $$\begin{array}{}& {\displaystyle }{f}_{\mathrm{hv}}=& {\displaystyle }\left(C_{\mathrm{v},N},,+,,C_{\mathrm{v},S}\right)\sqrt{C_{\mathrm{D}}}/k\left(R_i,>,\mathrm{0}\right),& {\displaystyle }=& {\displaystyle }\left(C_{\mathrm{v},N},,+,,C_{\mathrm{v},U}\right)\sqrt{C_{\mathrm{D}}}/k\left(R_i,\le ,\mathrm{0}\right),\end{array}$$ where the Richardson number $R_i$ is defined in Eq. (). The coefficients are defined as $$\begin{array}{}& {\displaystyle }{C}_{\mathrm{v},N}=& {\displaystyle }\ln \mathit{\{}\mathrm{1}+{\mathit{\varphi }}_z\left[\exp ,\left({\displaystyle \frac{k}{\sqrt{C_{\mathrm{DN}}}}}\right),-,\mathrm{1}\right]\mathit{\}},& {\displaystyle }{C}_{\mathrm{v},S}=& {\displaystyle }-{\mathit{\varphi }}_z\left({\displaystyle \frac{k}{\sqrt{C_{\mathrm{DN}}}}},-,{\displaystyle \frac{k}{\sqrt{C_{\mathrm{D}}}}}\right),& {\displaystyle }{C}_{\mathrm{v},U}=& {\displaystyle }-\ln \mathit{\{}\mathrm{1}+{\mathit{\varphi }}_z\left[\exp ,\left({\displaystyle \frac{k}{\sqrt{C_{\mathrm{DN}}}}},-,{\displaystyle \frac{k}{\sqrt{C_{\mathrm{D}}}}}\right),-,\mathrm{1}\right]\mathit{\}},\end{array}$$ where the drag coefficient $C_{\mathrm{D}}$ and the drag coefficient for neutral conditions $C_{\mathrm{DN}}$ are computed between the canopy air space and the free atmosphere above using the standard ISBA surface-layer transfer functions .

Longwave radiative flux expressions

The complete expression for the vegetation canopy net longwave radiation with an infinite number of reflections can be expressed as a series expansion e.g., as a function of the temperatures of the emitting surfaces ($T_{\mathrm{v}}$, $T_{\mathrm{g},\mathrm{1}}$, $T_{\mathrm{n},\mathrm{1}}$), their respective emissivities (${\mathit{ϵ}}_{\mathrm{v}}$, ${\mathit{ϵ}}_{\mathrm{g}}$ and ${\mathit{ϵ}}_{\mathrm{n}}$) and the canopy longwave absorption function, ${\mathit{\sigma }}_{\mathrm{LW}}$ (Eq. ). The MEB expressions are derived by explicitly expanding the series and assuming one reflection from each emitting source, which is a good approximation since emissivities are generally close to unity (fluxes from a single reflection are proportional to $\mathrm{1}-{\mathit{ϵ}}_x$ where $x$ represents g, v, or n, and $\mathit{ϵ}$ is close to unity for most natural surfaces).

Simple schematic for longwave radiation transfer for one reflection and up to three emitting surfaces (in addition to the down-welling atmospheric flux). Hollow arrows indicate fluxes after one reflection.

[Figure omitted. See PDF]

Snow is considered to be intercepted by the vegetation canopy and to accumulate on the ground below. The corresponding schematic of the radiative transfer is shown in Fig. . The canopy-intercepted snow is treated using a composite approach so that the canopy temperature $T_{\mathrm{v}}$ represents the effective temperature of the canopy-intercepted snow composite. The canopy emissivity is therefore simply defined as $${\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{v}}=\left(\mathrm{1},-,p_{\mathrm{nv}}\right){\mathit{ϵ}}_{\mathrm{v}}+p_{\mathrm{nv}}{\mathit{ϵ}}_{\mathrm{n}}.$$ In order to facilitate the use of a distinct multi-layer snow-process scheme, we split the fluxes between those interacting with the snowpack and the snow-free ground. The expressions for the snow-free surface are $$\begin{array}{}& {\displaystyle }& {\displaystyle }{A}_{\mathrm{g}}=\mathrm{LW}\downarrow \left(\mathrm{1},-,p_{\mathrm{ng}}\right),& {\displaystyle }& {\displaystyle }{B}_{\mathrm{g}}=A_{\mathrm{g}}{\mathit{\sigma }}_{\mathrm{LW}}\left(\mathrm{1},-,{\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{v}}\right),& {\displaystyle }& {\displaystyle }{C}_{\mathrm{g}}=A_{\mathrm{g}}\left(\mathrm{1},-,{\mathit{\sigma }}_{\mathrm{LW}}\right),& {\displaystyle }& {\displaystyle }{D}_{\mathrm{g}}=C_{\mathrm{g}}\left(\mathrm{1},-,{\mathit{ϵ}}_{\mathrm{g}}\right),& {\displaystyle }& {\displaystyle }{E}_{\mathrm{g}}=D_{\mathrm{g}}\left(\mathrm{1},-,{\mathit{\sigma }}_{\mathrm{LW}}^{\prime }\right),& {\displaystyle }& {\displaystyle }{F}_{\mathrm{g}}={\mathit{\sigma }}_{\mathrm{LW}}^{\prime }\mathit{\sigma }{\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{v}}{T}_{\mathrm{v}}^{\mathrm{4}}\left(\mathrm{1},-,p_{\mathrm{ng}}\right),& {\displaystyle }& {\displaystyle }{G}_{\mathrm{g}}=F_{\mathrm{g}}\left(\mathrm{1},-,{\mathit{ϵ}}_{\mathrm{g}}\right),& {\displaystyle }& {\displaystyle }{H}_{\mathrm{g}}=G_{\mathrm{g}}\left(\mathrm{1},-,{\mathit{\sigma }}_{\mathrm{LW}}^{\prime }\right),& {\displaystyle }& {\displaystyle }{I}_{\mathrm{g}}=\mathit{\sigma }{\mathit{ϵ}}_{\mathrm{g}}{T}_{\mathrm{g}}^{\mathrm{4}}\left(\mathrm{1},-,p_{\mathrm{ng}}\right),& {\displaystyle }& {\displaystyle }{J}_{\mathrm{g}}=I_{\mathrm{g}}{\mathit{\sigma }}_{\mathrm{LW}}^{\prime }\left(\mathrm{1},-,{\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{v}}\right)\left(\mathrm{1},-,p_{\mathrm{ng}}^{\prime }\right),& {\displaystyle }& {\displaystyle }{K}_{\mathrm{g}}=I_{\mathrm{g}}{\mathit{\sigma }}_{\mathrm{LW}}^{\prime }\left(\mathrm{1},-,{\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_{\mathrm{v}}\right)p_{\mathrm{ng}}^{\prime },& {\displaystyle }& {\displaystyle }{L}_{\mathrm{g}}=I_{\mathrm{g}}\left(\mathrm{1},-,{\mathit{\sigma }}_{\mathrm{LW}}^{\prime }\right),& {\displaystyle }& {\displaystyle }{p}_{\mathrm{ng}}^{\prime }=p_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right),\end{array}$$ and the equations for the snow-covered understory fraction are $$\begin{array}{}& {\displaystyle }& {\displaystyle }{A}_{\mathrm{n}}=\mathrm{LW}\downarrow p_{\mathrm{ng}},& {\displaystyle }& {\displaystyle }{B}_{\mathrm{n}}=A_{\mathrm{n}}{\mathit{\sigma }}_{\mathrm{f}\mathrm{LW}}\left(\mathrm{1},-,{\mathit{ϵ}}_{\mathrm{v}}\right),& {\displaystyle }& {\displaystyle }{C}_{\mathrm{n}}=A_{\mathrm{n}}\left(\mathrm{1},-,{\mathit{\sigma }}_{\mathrm{f}\mathrm{LW}}\right),& {\displaystyle }& {\displaystyle }{D}_{\mathrm{n}}=C_{\mathrm{n}}\left(\mathrm{1},-,{\mathit{ϵ}}_{\mathrm{n}}\right),& {\displaystyle }& {\displaystyle }{E}_{\mathrm{n}}=D_{\mathrm{n}}\left(\mathrm{1},-,{\mathit{\sigma }}_{\mathrm{LW}}^{\prime }\right),& {\displaystyle }& {\displaystyle }{F}_{\mathrm{n}}={\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{f}\mathrm{LW}}\mathit{\sigma }{\mathit{ϵ}}_{\mathrm{v}}{T}_{\mathrm{v}}^{\mathrm{4}}{p}_{\mathrm{ng}},& {\displaystyle }& {\displaystyle }{G}_{\mathrm{n}}=F_{\mathrm{n}}\left(\mathrm{1},-,{\mathit{ϵ}}_{\mathrm{n}}\right),& {\displaystyle }& {\displaystyle }{H}_{\mathrm{n}}=G_{\mathrm{n}}\left(\mathrm{1},-,{\mathit{\sigma }}_{\mathrm{LW}}^{\prime }\right),& {\displaystyle }& {\displaystyle }{I}_{\mathrm{n}}=\mathit{\sigma }{\mathit{ϵ}}_{\mathrm{n}}{T}_{\mathrm{n}}^{\mathrm{4}}{p}_{\mathrm{ng}},& {\displaystyle }& {\displaystyle }{J}_{\mathrm{n}}=I_{\mathrm{n}}{\mathit{\sigma }}_{\mathrm{LW}}^{\prime }\left(\mathrm{1},-,{\mathit{ϵ}}_{\mathrm{v}}\right)\left(\mathrm{1},-,p_{\mathrm{ng}}^{\prime \prime }\right),& {\displaystyle }& {\displaystyle }{K}_{\mathrm{n}}=I_{\mathrm{n}}{\mathit{\sigma }}_{\mathrm{LW}}^{\prime }\left(\mathrm{1},-,{\mathit{ϵ}}_{\mathrm{v}}\right)p_{\mathrm{ng}}^{\prime \prime },& {\displaystyle }& {\displaystyle }{L}_{\mathrm{n}}=I_{\mathrm{n}}\left(\mathrm{1},-,{\mathit{\sigma }}_{\mathrm{LW}}^{\prime }\right),& {\displaystyle }& {\displaystyle }{p}_{\mathrm{ng}}^{\prime \prime }=p_{\mathrm{ng}}+p_{\mathrm{n}\mathit{\alpha }}\left(\mathrm{1},-,p_{\mathrm{ng}}\right),\end{array}$$ where the different terms are indicated in Fig. . In MEB, the ground-based snowpack depth can increase to the point that it buries the canopy; thus, for both the snow-covered and snow-free understory fractions a modified snow fraction is defined as $${\mathit{\sigma }}_{\mathrm{LW}}^{\prime }=\left(\mathrm{1},-,p_{\mathrm{ng}}^{\prime }\right){\mathit{\sigma }}_{\mathrm{LW}}+p_{\mathrm{ng}}^{\prime }{\mathit{\sigma }}_{\mathrm{f}\mathrm{LW}}.$$ The factor, ${\mathit{\sigma }}_{\mathrm{f}\mathrm{LW}}$, over the understory snow-covered fraction of the grid box is modeled quite simply from Eq. (). The net longwave radiation for the understory, snowpack, and vegetation canopy are therefore defined, respectively, as $$\begin{array}{cc}& {\displaystyle }{\mathrm{LW}}_{\mathrm{net}\mathrm{g}}=& {\displaystyle }{C}_{\mathrm{g}}+F_{\mathrm{g}}+J_{\mathrm{g}}+J_{\mathrm{n}}-D_{\mathrm{g}}-G_{\mathrm{g}}-I_{\mathrm{g}},& {\displaystyle }{\mathrm{LW}}_{\mathrm{net}\mathrm{n}}=& {\displaystyle }{C}_{\mathrm{n}}+F_{\mathrm{n}}+K_{\mathrm{n}}+K_{\mathrm{g}}-D_{\mathrm{n}}-G_{\mathrm{n}}-I_{\mathrm{n}},{\displaystyle }{\mathrm{LW}}_{\mathrm{net}\mathrm{v}}=& {\displaystyle }{A}_{\mathrm{g}}+D_{\mathrm{g}}+G_{\mathrm{g}}+I_{\mathrm{g}}+A_{\mathrm{n}}+D_{\mathrm{n}}+G_{\mathrm{n}}\\ {\displaystyle }& {\displaystyle }+I_{\mathrm{n}}-B_{\mathrm{g}}-C_{\mathrm{g}}-E_{\mathrm{g}}-H_{\mathrm{g}}-\mathrm{\hspace{0.17em}2}{F}_{\mathrm{g}}\\ {\displaystyle }& {\displaystyle }-J_{\mathrm{g}}-L_{\mathrm{g}}-K_{\mathrm{g}}-B_{\mathrm{n}}-C_{\mathrm{n}}-E_{\mathrm{n}}\\ & {\displaystyle }& {\displaystyle }-H_{\mathrm{n}}-\mathrm{\hspace{0.17em}2}{F}_{\mathrm{n}}-J_{\mathrm{n}}-L_{\mathrm{n}}-K_{\mathrm{n}},\end{array}$$ where the upwelling longwave radiation is computed from $$\mathrm{LW}\uparrow =\mathrm{LW}\downarrow -{\mathrm{LW}}_{\mathrm{net}\mathrm{g}}-{\mathrm{LW}}_{\mathrm{net}\mathrm{n}}-{\mathrm{LW}}_{\mathrm{net}\mathrm{v}}.$$ The inclusion of the snow-buried canopy fraction in Eqs. () and () causes all of the vegetation transmission and below canopy fluxes to vanish as $p_{\mathrm{ng}}$ and $p_{\mathrm{n}\mathit{\alpha }}\to \mathrm{0}$ so that the only longwave radiative exchanges occur between the atmosphere and the snowpack in this limit.

Net longwave radiation flux derivatives

The first-order derivatives of the net longwave radiation terms are needed in order to solve the system of linearized surface energy budget equations (Eqs. –). The Taylor series expansion (neglecting higher-order terms) is expressed as $$\begin{array}{cc}{\displaystyle }{\mathrm{LW}}_{\mathrm{net}\mathrm{i}}^{+}=& {\displaystyle }{\mathrm{LW}}_{\mathrm{net}\mathrm{i}}+\sum _{j=\mathrm{1}}^{N_{\mathrm{seb}}}{\displaystyle \frac{\partial L_{\mathrm{net}\mathrm{i}}}{\partial T_j}}\left(T_j^{+},-,T_j\right)\\ & {\displaystyle }& {\displaystyle }\left(i,=,\mathrm{1},,,N_{\mathrm{seb}}\right),\end{array}$$ where $N_{\mathrm{seb}}$ represents the number of surface energy budgets, and $i$ and $j$ represent the indexes for each energy budget. The superscript $+$ represents the variable at time $t+\mathrm{\Delta }t$, while by default, no superscript represents the value at time $t$. Equation () therefore results in a $N_{\mathrm{seb}}\times N_{\mathrm{seb}}$ Jacobian matrix ($\mathrm{3}\times \mathrm{3}$ for MEB). The matrix coefficients are expressed as $$\begin{array}{cc}{\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{v}}}{\partial T_{\mathrm{v}}}}=& {\displaystyle \frac{\partial G_{\mathrm{g}}}{\partial T_{\mathrm{v}}}}-{\displaystyle \frac{\partial H_{\mathrm{g}}}{\partial T_{\mathrm{v}}}}-\mathrm{\hspace{0.17em}2}{\displaystyle \frac{\partial F_{\mathrm{g}}}{\partial T_{\mathrm{v}}}}+{\displaystyle \frac{\partial G_{\mathrm{n}}}{\partial T_{\mathrm{v}}}}-{\displaystyle \frac{\partial H_{\mathrm{n}}}{\partial T_{\mathrm{v}}}}\\ & {\displaystyle }& {\displaystyle }-\mathrm{\hspace{0.17em}2}{\displaystyle \frac{\partial F_{\mathrm{n}}}{\partial T_{\mathrm{v}}}},& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{v}}}{\partial T_{\mathrm{g}}}}=& {\displaystyle \frac{\partial I_{\mathrm{g}}}{\partial T_{\mathrm{g}}}}-{\displaystyle \frac{\partial J_{\mathrm{g}}}{\partial T_{\mathrm{g}}}}-{\displaystyle \frac{\partial K_{\mathrm{g}}}{\partial T_{\mathrm{g}}}}-{\displaystyle \frac{\partial L_{\mathrm{g}}}{\partial T_{\mathrm{g}}}},& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{v}}}{\partial T_{\mathrm{n}}}}=& {\displaystyle \frac{\partial I_{\mathrm{n}}}{\partial T_{\mathrm{n}}}}-{\displaystyle \frac{\partial J_{\mathrm{n}}}{\partial T_{\mathrm{n}}}}-{\displaystyle \frac{\partial K_{\mathrm{n}}}{\partial T_{\mathrm{n}}}}-{\displaystyle \frac{\partial L_{\mathrm{n}}}{\partial T_{\mathrm{n}}}},& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{g}}}{\partial T_{\mathrm{v}}}}=& {\displaystyle \frac{\partial F_{\mathrm{g}}}{\partial T_{\mathrm{v}}}}-{\displaystyle \frac{\partial G_{\mathrm{g}}}{\partial T_{\mathrm{v}}}},& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{g}}}{\partial T_{\mathrm{g}}}}=& {\displaystyle \frac{\partial J_{\mathrm{g}}}{\partial T_{\mathrm{g}}}}-{\displaystyle \frac{\partial I_{\mathrm{g}}}{\partial T_{\mathrm{g}}}},& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{g}}}{\partial T_{\mathrm{n}}}}=& {\displaystyle \frac{\partial J_{\mathrm{n}}}{\partial T_{\mathrm{n}}}},& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{n}}}{\partial T_{\mathrm{v}}}}=& {\displaystyle \frac{\partial F_{\mathrm{n}}}{\partial T_{\mathrm{v}}}}-{\displaystyle \frac{\partial G_{\mathrm{n}}}{\partial T_{\mathrm{v}}}},& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{n}}}{\partial T_{\mathrm{g}}}}=& {\displaystyle \frac{\partial K_{\mathrm{g}}}{\partial T_{\mathrm{g}}}},& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{n}}}{\partial T_{\mathrm{n}}}}=& {\displaystyle \frac{\partial J_{\mathrm{n}}}{\partial T_{\mathrm{n}}}}-{\displaystyle \frac{\partial I_{\mathrm{n}}}{\partial T_{\mathrm{n}}}}.\end{array}$$ Using Eq. () to evaluate the derivatives we have $$\begin{array}{}& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{v}}}{\partial T_{\mathrm{v}}}}=& {\displaystyle \frac{\mathrm{4}}{T_{\mathrm{v}}}}\left(G_{\mathrm{g}},-,H_{\mathrm{g}},-,\mathrm{2},F_{\mathrm{g}},+,G_{\mathrm{n}},-,H_{\mathrm{n}},-,\mathrm{2},F_{\mathrm{n}}\right),& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{v}}}{\partial T_{\mathrm{g}}}}=& {\displaystyle \frac{\mathrm{4}}{T_{\mathrm{g}}}}\left(I_{\mathrm{g}},-,J_{\mathrm{g}},-,K_{\mathrm{g}},-,L_{\mathrm{g}}\right),& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{v}}}{\partial T_{\mathrm{n}}}}=& {\displaystyle \frac{\mathrm{4}}{T_{\mathrm{n}}}}\left(I_{\mathrm{n}},-,J_{\mathrm{n}},-,K_{\mathrm{n}},-,L_{\mathrm{n}}\right),& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{g}}}{\partial T_{\mathrm{v}}}}=& {\displaystyle \frac{\mathrm{4}}{T_{\mathrm{v}}}}\left(F_{\mathrm{g}},-,G_{\mathrm{g}}\right),& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{g}}}{\partial T_{\mathrm{g}}}}=& {\displaystyle \frac{\mathrm{4}}{T_{\mathrm{g}}}}\left(J_{\mathrm{g}},-,I_{\mathrm{g}}\right),& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{g}}}{\partial T_{\mathrm{n}}}}=& {\displaystyle \frac{\mathrm{4}}{T_{\mathrm{n}}}}{J}_{\mathrm{n}},& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{n}}}{\partial T_{\mathrm{v}}}}=& {\displaystyle \frac{\mathrm{4}}{T_{\mathrm{v}}}}\left(F_{\mathrm{n}},-,G_{\mathrm{n}}\right),& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{n}}}{\partial T_{\mathrm{g}}}}=& {\displaystyle \frac{\mathrm{4}}{T_{\mathrm{g}}}}{K}_{\mathrm{g}},& {\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{n}}}{\partial T_{\mathrm{n}}}}=& {\displaystyle \frac{\mathrm{4}}{T_{\mathrm{n}}}}\left(J_{\mathrm{n}},-,I_{\mathrm{n}}\right),\end{array}$$ and therefore from a coding perspective, the computation of the derivatives is trivial (using already computed quantities).

Halstead coefficient maximum

A maximum Halstead coefficient is imposed by estimating which value of ${\mathit{\delta }}_{\mathrm{v}}$ that is needed to just evaporate any existing intercepted water $W_{\mathrm{rv}}$ given the conditions at the beginning of the time step. Assuming that phase changes are small, and neglecting canopy drip and any condensation from transpiration, the time-differenced prognostic equation for intercepted water on canopy vegetation (Eq. ) can be approximated as $${\displaystyle \frac{{W_{\mathrm{rv}}}^{+}-W_{\mathrm{rv}}}{\mathrm{\Delta }t}}=(\mathrm{1}-{\mathit{\chi }}_{\mathrm{v}})(\mathrm{1}-p_{\mathrm{ng}}{p}_{\mathit{\alpha }\mathrm{n}})P_{\mathrm{r}}-E_{\mathrm{r}}.$$ Assuming that all existing water evaporates in one time step (i.e., $W_{\mathrm{rv}}^{+}=\mathrm{0}$), and substituting the full expression for $E_{\mathrm{r}}$ (Eq. ) into Eq. (), the maximum value of ${\mathit{\delta }}_{\mathrm{v}}$ can be determined as

$$\begin{array}{cc}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}_{\mathrm{v},\max }={\displaystyle }& {\displaystyle \frac{\left[\left(\mathrm{1},-,{\mathit{\chi }}_{\mathrm{v}}\right),\left(\mathrm{1},-,p_{\mathrm{ng}},p_{\mathit{\alpha }\mathrm{n}}\right),P_{\mathrm{r}},,+,,\left(W_{\mathrm{rv}},/,\mathrm{\Delta },t\right)\right]\left(L,/,L_{\mathrm{v}}\right)}{{\mathit{\rho }}_{\mathrm{a}}\left(\mathrm{1},-,p_{\mathrm{nv}}\right)k_{\mathrm{v}}\mathit{\{}\left[p_{\mathrm{ng}},\left(\mathrm{1},-,p_{\mathit{\alpha }\mathrm{n}}\right),/,R_{\mathrm{avn}-\mathrm{c}}\right]+\left[\left(\mathrm{1},-,p_{\mathrm{ng}}\right),/,R_{\mathrm{avg}-\mathrm{c}}\right]\mathit{\}}\left(q_{\mathrm{sat}\mathrm{v}},,-,,q_{\mathrm{c}}\right)}}.\end{array}$$ Equation () is an approximation since all of the variables on the RHS use conditions from the start of the time step; however, this method has proven to greatly reduce the risk for occasional numerical artifacts (jumps) and the associated need for mass corrections (if net losses in mass exceed the updated test value for interception storage).

Energy and mass conservation

Energy conservation

The soil and snowpack prognostic temperature equations can be written in flux form for $k=\mathrm{1},N_{\mathrm{g}}$ soil layers and $k=\mathrm{1},N_{\mathrm{n}}$ snow layers as $$\begin{array}{cc}& {\displaystyle }{\mathcal{C}}_{\mathrm{g},k}{\displaystyle \frac{\partial T_{\mathrm{g},k}}{\partial t}}& {\displaystyle }=G_{\mathrm{g},k-\mathrm{1}}-G_{\mathrm{g},k}+L_{\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{g},k},{\displaystyle }{\mathcal{C}}_{\mathrm{n},k}{\displaystyle \frac{\partial T_{\mathrm{n},k}}{\partial t}}& {\displaystyle }=G_{\mathrm{n},k-\mathrm{1}}-G_{\mathrm{n},k}+L_{\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{n},k}+{\mathit{\xi }}_{\mathrm{n},k-\mathrm{1}}\\ & {\displaystyle }& {\displaystyle }-{\mathit{\xi }}_{\mathrm{n},k}+{\mathrm{SW}}_{\mathrm{net},\mathrm{n}}\left({\mathit{\tau }}_{\mathrm{n},k-\mathrm{1}},-,{\mathit{\tau }}_{\mathrm{n},k}\right).\end{array}$$ The total energy balance of the vegetation canopy–soil–snowpack system is conserved at each time step $\mathrm{\Delta }t$ and can be obtained by summing the discrete time forms of Eqs. (), (), and () for the vegetation and all soil and snow layers, respectively, yielding $${\displaystyle }\begin{array}{l}{\mathcal{C}}_{\mathrm{v}}\mathrm{\Delta }{T}_{\mathrm{v}}+\sum _{k=\mathrm{1}}^{N_{\mathrm{g}}}{\mathcal{C}}_{\mathrm{g},k}\mathrm{\Delta }{T}_{\mathrm{g},k}+p_{\mathrm{ng}}\sum _{k=\mathrm{1}}^{N_{\mathrm{n}}}{\mathcal{C}}_{\mathrm{n},k}\mathrm{\Delta }{T}_{\mathrm{n},k}=\\ \mathrm{\Delta }t[\left(\mathrm{1},-,p_{\mathrm{ng}}\right)G_{\mathrm{g},\mathrm{0}}++p_{\mathrm{ng}}\left(G_{\mathrm{gn}},,+,,{\mathit{\tau }}_{\mathrm{n},N_{\mathrm{n}}},{\mathrm{SW}}_{\mathrm{net},\mathrm{n}},,+,,G_{\mathrm{n},\mathrm{0}}\right),\\ +R_{\mathrm{n}\mathrm{v}}-H_{\mathrm{v}}-{\mathrm{LE}}_{\mathrm{v}}+L_{\mathrm{f}}\left({\mathrm{\Phi }}_{\mathrm{v}},,+,,\sum _{k=\mathrm{1}}^{N_{\mathrm{g}}},{\mathrm{\Phi }}_{\mathrm{g},k}\right.\\ \left..,+,,p_{\text{ng}},\sum _{k=\mathrm{1}}^{N_{\mathrm{n}}},{\mathrm{\Phi }}_{\mathrm{n},k}\right)],\end{array}$$ where $\mathrm{\Delta }{T}_x=T_x(t+\mathrm{\Delta }t)-T_x\left(t\right)$. Note that Eq. () must first be multiplied by $p_{\mathrm{ng}}$ in order to make it patch or grid cell relative when it is combined with the soil and vegetation budget equations. The surface boundary conditions for Eqs. () and () are, respectively, $$\begin{array}{}& {\displaystyle }{G}_{\mathrm{g},\mathrm{0}}& {\displaystyle }=R_{\mathrm{n}\mathrm{g}}-H_{\mathrm{g}}-{\mathrm{LE}}_{\mathrm{g}},& {\displaystyle }{G}_{\mathrm{n},\mathrm{0}}& {\displaystyle }=R_{\mathrm{n}\mathrm{n}}-H_{\mathrm{n}}-{\mathrm{LE}}_{\mathrm{n}},& {\displaystyle }{\mathit{\tau }}_{\mathrm{n},\mathrm{0}}& {\displaystyle }=\mathrm{1},& {\displaystyle }{\mathit{\xi }}_{\mathrm{n},\mathrm{0}}& {\displaystyle }=0.\end{array}$$ Equation () signifies that the net shortwave radiation at the surface enters the snowpack, and Eq. () represents the fact that energy changes owing to the time-evolving snow grid can only arise in the surface layer owing to exchanges with the sub-surface layer. Snowfall is assumed to have the same temperature as the snowpack; thus, a corresponding cooling/heating term does not appear in Eq. (), although the corresponding mass increase must appear in the snow water budget equation (see Sect. ).

The lower boundary conditions for Eqs. () and () are, respectively, $$\begin{array}{}& {\displaystyle }{G}_{\mathrm{g},N_{\mathrm{g}}}& {\displaystyle }=\mathrm{0},& {\displaystyle }{\mathit{\xi }}_{\mathrm{n},N_{\mathrm{n}}}& {\displaystyle }=0.\end{array}$$

The appearance of the same discrete form for $\mathrm{\Phi }$ in both the energy and mass budget equations ensures enthalpy conservation. Owing to Eqs. () and (), the total effective heating of the snowpack owing to grid adjustments is $$\underset{\mathrm{0}}{\overset{D_{N_{\mathrm{n}}}}{\int }}{\mathit{\xi }}_{\mathrm{n}}d{D}_{\mathrm{n}}=\mathrm{\hspace{0.17em}\hspace{0.17em}0},$$ where $D_{N_{\mathrm{n}}}$ represents the total snow depth. Thus, this term only represents a contribution from contiguous snow layers, not from a source external to the snowpack. The energy storage of the snow–soil–vegetation system is balanced by the net surface radiative and turbulent fluxes and internal phase changes (solid and liquid phases of water substance).

Mass conservation

The soil and snowpack prognostic mass equations can be written in flux form for $k=\mathrm{2},N_{\mathrm{g}\mathrm{w}}$ soil layers and $k=\mathrm{1},N_{\mathrm{n}}$ snow layers as $$\begin{array}{cc}{\displaystyle }& {\displaystyle \frac{\partial W_{\mathrm{n},k}}{\partial t}}=F_{\mathrm{nl},k-\mathrm{1}}-F_{\mathrm{nl},k}-{\mathrm{\Phi }}_{\mathrm{n},k}\\ & {\displaystyle }& {\displaystyle }+{\mathit{\xi }}_{\mathrm{nl},k}-{\mathit{\xi }}_{\mathrm{nl},k-\mathrm{1}}\left(k,=,\mathrm{2},,,N_{\mathrm{n}}\right),{\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\mathrm{w}}\mathrm{\Delta }{z}_{\mathrm{g},\mathrm{1}}{\displaystyle \frac{\partial w_{\mathrm{g},k}}{\partial t}}=F_{\mathrm{g},k-\mathrm{1}}-F_{\mathrm{g},k}-{\mathrm{\Phi }}_{\mathrm{g},k}\\ & {\displaystyle }& {\displaystyle }-{\mathcal{F}}_{\mathrm{2},k}\max \left(\mathrm{0},,,,E_{\mathrm{tr}}\right)\left(k,=,\mathrm{2},,,N_{\mathrm{g}\mathrm{w}}\right),& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\mathrm{w}}\mathrm{\Delta }{z}_{\mathrm{g},\mathrm{1}}{\displaystyle \frac{\partial w_{\mathrm{gf},k}}{\partial t}}={\mathrm{\Phi }}_{\mathrm{g},k}\left(k,=,\mathrm{2},,,N_{\mathrm{g}\mathrm{w}}\right).\end{array}$$ The total grid box water budget at each time step is obtained by summing the budget equations for the surface layers (Eqs. –) together with those for the sub-surface layers (Eqs. –) to have $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathrm{\Delta }{W}_{\mathrm{r}}+\mathrm{\Delta }{W}_{\mathrm{r}\mathrm{n}}+p_{\mathrm{ng}}\sum _{k=\mathrm{1}}^{N_{\mathrm{n}}}\mathrm{\Delta }{W}_{\mathrm{n},k}+{\mathit{\rho }}_{\mathrm{w}}\sum _{k=\mathrm{1}}^{N_{\mathrm{gw}}}\mathrm{\Delta }{z}_{g,k}\left(w_{\mathrm{g}k},+,w_{\mathrm{gf}k}\right)\\ {\displaystyle }& {\displaystyle }=\mathrm{\Delta }t[P_{\mathrm{r}}+P_{\mathrm{s}}-R_{\mathrm{0}}-F_{\mathrm{g},N_{\mathrm{gw}}}-\left(\mathrm{1},-,p_{\mathrm{ng}}\right)E_{\mathrm{g}}-E_{\mathrm{v}}\\ & {\displaystyle }& {\displaystyle }-p_{\mathrm{ng}}{E}_{\mathrm{n}}-{\mathrm{\Phi }}_{\mathrm{v}}-\sum _{k=\mathrm{1}}^{N_{\mathrm{g}}}{\mathrm{\Phi }}_{\mathrm{g},k}-p_{\mathrm{ng}}\sum _{k=\mathrm{1}}^{N_{\mathrm{n}}}{\mathrm{\Phi }}_{\mathrm{n},k}],\end{array}$$ where Eq. () has been multiplied by $p_{\mathrm{ng}}$ to make it patch or grid box relative (as was done for energy conservation in Sect. ). $R_{\mathrm{0}}$ can simply be a diagnostic or coupled with a river routing scheme . The soil water lower boundary condition $F_{\mathrm{g},N_{\mathrm{gw}}}$ represents the base flow or drainage leaving the lowest hydrological layer, which can then be transferred as input to a river routing scheme (see references above) or to a ground water scheme. In such instances, it can be negative if an option to permit a ground water inflow is activated . The soil liquid water and equivalent frozen water equivalent volumetric water content extend down to layer $N_{\mathrm{g}\mathrm{w}}$, where $N_{\mathrm{g}\mathrm{w}}\le N_{\mathrm{g}}$. Note that the vertical soil water transfer or evolution is not computed below $z_{\mathrm{g}}\left(k,=,N_{\mathrm{gw}}\right)$, whereas heat transfer can be. In order to compute the thermal properties for deep soil temperature (thermal conductivity and heat capacity for example), soil moisture estimates are needed: values from the soil are extrapolated downward assuming hydrostatic equilibrium A detailed description of the soil model is given by and .

Note that Eq. () is snow relative; therefore, this equation must be multiplied by the ground-based snow fraction $p_{\mathrm{ng}}$ to be grid box relative for coupling with the soil and vegetation water storage terms. The lower boundary condition for liquid water flow $F_{\mathrm{nl},N_{\mathrm{n}}}$ is defined as the liquid water exceeding the lowest maximum snow-layer liquid water-holding capacity. ${\mathit{\xi }}_{\mathrm{nl}}$ represents the internal mass changes of a snowpack layer when the vertical grid is reset. When integrated over the entire snowpack depth, this term vanishes (analogous to Eq. () for the snowpack temperature equation). See and for details on the snow model processes.

The equations describing flooding are not described in detail here as this parameterization is independent of MEB, and it is described in detail by . The coupling of MEB with the interactive flooding scheme will be the subject of a future paper.

Implicit numerical coupling with the atmosphere

The land–atmosphere coupling is accomplished through the atmospheric model vertical diffusion (heat, mass, momentum, chemical species, aerosols, etc.) and radiative schemes. Owing to the potential for relatively large diffusivity, especially in the lower atmosphere near the surface, fairly strict time step constraints must be applied. In this section, a fully implicit time scheme (with an option for explicit coupling) is described. There are two reasons for using this approach; (i) an implicit coupling is more numerically stable not only for time steps typical of GCM applications but also for some NWP models, and (ii) the methodology permits code modularity in that the land surface model routines can be independent of the atmospheric model code and they can be called using a standard interface, which is the philosophy of SURFEX . The coupling follows the methodology first proposed by , which was further generalized by .

The atmospheric turbulence scheme is generally expressed as a second-order diffusion equation in the vertical (which is assumed herein) and it is discretized using the backward difference time scheme. Note that a semi-implicit scheme, such as the Crank–Nicolson , could also be used within this framework. Thus, the equations can be cast as a tri-diagonal matrix. Assuming a fixed for zero (the general case) upper boundary condition at the top of the atmosphere, the diffusion equations for the generic variable $\mathit{\varphi }$ can be cast as a linear function of the variable in the layer below as $${\mathit{\varphi }}_k^{+}=B_{\mathit{\varphi },k}+A_{\mathit{\varphi },k}{\mathit{\varphi }}_{k+\mathrm{1}}^{+}\left(k,=,\mathrm{1},,,N_{\mathrm{a}},-,\mathrm{1}\right),$$ where $N_{\mathrm{a}}$ represents the number of atmospheric model layers, $k=\mathrm{1}$ represents the uppermost layer with $k$ increasing with decreasing height above the surface, and the superscript $+$ indicates the value of $\mathit{\varphi }$ at time $t+\mathrm{\Delta }t$ (at the end of the time step). The coefficients $A_{\mathit{\varphi },k}$ and $B_{\mathit{\varphi },k}$ are computed in a downward sweep within the turbulence scheme and thus consist in atmospheric prognostic variables, diffusivity, heat capacities, and additional source terms from layer $k$ and above evaluated at time level $t$ . As shown by , the equation for the lowest atmospheric model layer can be expressed using a flux lower boundary condition as $${\mathit{\varphi }}_{N_{\mathrm{a}}}^{+}=B_{\mathit{\varphi },N_{\mathrm{a}}}+A_{\mathit{\varphi },N_{\mathrm{a}}}{F}_{\mathit{\varphi },N_{\mathrm{a}}+\mathrm{1}}^{+},$$ where $F_{\mathit{\varphi },N_{\mathrm{a}}+\mathrm{1}}^{+}$ is the implicit surface flux from one or multiple surface energy budgets. Technically, only the $B_{\mathit{\varphi },N_{\mathrm{a}}}$ and $A_{\mathit{\varphi },N_{\mathrm{a}}}$ coefficients are needed by the LSM in order to compute the updated land surface fluxes and temperatures, which are fully implicitly coupled with the atmosphere. Once $F_{\mathit{\varphi },N_{\mathrm{a}}+\mathrm{1}}^{+}$ has been computed by the LSM, it can be returned to the atmospheric turbulence scheme, which can then solve for ${\mathit{\varphi }}_k^{+}$ from $k=N_{\mathrm{a}}$ to $k=\mathrm{1}$ (i.e., the upward sweep). For explicit land–atmosphere coupling or offline land-only applications, the coupling coefficients can be set to $A_{\mathit{\varphi },N_{\mathrm{a}}}=\mathrm{0}$ and $B_{\mathit{\varphi },N_{\mathrm{a}}}={\mathit{\varphi }}_{N_{\mathrm{a}}}$ in the driving code.

Numerical solution of the surface energy budgets

Discretization of surface energy budgets

The surface energy budget equations (Eqs. –) are integrated in time using the implicit backward difference scheme. They can be written in discretized form as $$\begin{array}{cc}{\displaystyle }{\mathcal{C}}_{\mathrm{v}}& {\displaystyle \frac{\left(T_{\mathrm{v}}^{+},-,T_{\mathrm{v}}\right)}{\mathrm{\Delta }t}}={\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{v}}}{\partial T_{\mathrm{v}}}}\\ {\displaystyle }& {\displaystyle }\left(T_{\mathrm{v}}^{+},-,T_{\mathrm{v}}\right)+{\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{v}}}{\partial T_{\mathrm{g},\mathrm{1}}}}\left(T_{\mathrm{g},\mathrm{1}}^{+},-,T_{\mathrm{g},\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{v}}}{\partial T_{\mathrm{n},\mathrm{1}}}}\left(T_{\mathrm{n},\mathrm{1}}^{+},-,T_{\mathrm{n},\mathrm{1}}\right)+{\mathrm{SW}}_{\mathrm{net}\mathrm{v}}+{\mathrm{LW}}_{\mathrm{net}\mathrm{v}}\\ {\displaystyle }& {\displaystyle }+{\mathit{\phi }}_{\mathrm{v}}\left({\mathcal{A}}_{\mathrm{v}},,T_{\mathrm{v}}^{+},-,{\mathcal{A}}_{\mathrm{c}},,T_{\mathrm{c}}^{+}\right)\\ & {\displaystyle }& {\displaystyle }+h_{\mathrm{sv}}{\mathit{\phi }}_{\mathrm{v}}L\left[q_{\mathrm{sat}\mathrm{v}},,+,,{\displaystyle \frac{\partial q_{\mathrm{sat}\mathrm{v}}}{\partial T_{\mathrm{v}}}},\left(T_{\mathrm{v}}^{+},-,T_{\mathrm{v}}\right),,-,,q_{\mathrm{c}}^{+}\right],{\displaystyle }& {\displaystyle }{\mathcal{C}}_{\mathrm{g},\mathrm{1}}{\displaystyle \frac{\left(T_{\mathrm{g},\mathrm{1}}^{+},-,T_{\mathrm{g},\mathrm{1}}\right)}{\mathrm{\Delta }t}}=[{\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{g}}}{\partial T_{\mathrm{v}}}}\\ {\displaystyle }& {\displaystyle }\left(T_{\mathrm{v}}^{+},-,T_{\mathrm{v}}\right)+{\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{g}}}{\partial T_{\mathrm{g},\mathrm{1}}}}\left(T_{\mathrm{g},\mathrm{1}}^{+},-,T_{\mathrm{g},\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{g}}}{\partial T_{\mathrm{n},\mathrm{1}}}}\left(T_{\mathrm{n},\mathrm{1}}^{+},-,T_{\mathrm{n},\mathrm{1}}\right)+{\mathrm{SW}}_{\mathrm{net}\mathrm{g}}+{\mathrm{LW}}_{\mathrm{net}\mathrm{g}}\\ {\displaystyle }& {\displaystyle }+{\mathit{\phi }}_{\mathrm{g}}\left({\mathcal{A}}_{\mathrm{g}},,T_{\mathrm{g}}^{+},-,{\mathcal{A}}_{\mathrm{c}},,T_{\mathrm{c}}^{+}\right)\\ {\displaystyle }& {\displaystyle }+{\mathit{\phi }}_{\mathrm{g}}L\mathit{\{}{h}_{\mathrm{sg}}\left[q_{\mathrm{sat}\mathrm{g}},,+,,{\displaystyle \frac{\partial q_{\mathrm{sat}\mathrm{g}}}{\partial T_{\mathrm{g}}}},\left(T_{\mathrm{g}}^{+},-,T_{\mathrm{g}}\right)\right]-h_{\mathrm{a}}{q}_{\mathrm{c}}^{+}\mathit{\}}\\ {\displaystyle }& {\displaystyle }]\left(\mathrm{1},-,p_{\mathrm{ng}}\right)+p_{\mathrm{ng}}{\mathrm{\Lambda }}_{\mathrm{g},\mathrm{n}}\left(T_{\mathrm{n},N_{\mathrm{n}}}^{\ast },-,T_{\mathrm{g},\mathrm{1}}^{+}\right),\\ & {\displaystyle }& {\displaystyle }-{\mathrm{\Lambda }}_{\mathrm{g},\mathrm{1}}\left(T_{\mathrm{g},\mathrm{1}}^{+},-,T_{\mathrm{g},\mathrm{2}}^{+}\right),{\displaystyle }& {\displaystyle }{p}_{\mathrm{ng}}{\mathcal{C}}_{\mathrm{n},\mathrm{1}}{\displaystyle \frac{\left(T_{\mathrm{n},\mathrm{1}}^{+},-,T_{\mathrm{n},\mathrm{1}}\right)}{\mathrm{\Delta }t}}=\mathit{\{}{\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{n}}}{\partial T_{\mathrm{v}}}}\left(T_{\mathrm{v}}^{+},-,T_{\mathrm{v}}\right)\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{n}}}{\partial T_{\mathrm{g},\mathrm{1}}}}\left(T_{\mathrm{g},\mathrm{1}}^{+},-,T_{\mathrm{g},\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial {\mathrm{LW}}_{\mathrm{net}\mathrm{n}}}{\partial T_{\mathrm{n},\mathrm{1}}}}\left(T_{\mathrm{n},\mathrm{1}}^{+},-,T_{\mathrm{n},\mathrm{1}}\right)+{\mathrm{SW}}_{\mathrm{net}\mathrm{n}}+{\mathrm{LW}}_{\mathrm{net}\mathrm{n}}\\ {\displaystyle }& {\displaystyle }+\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right){\mathit{\phi }}_{\mathrm{n}-\mathrm{c}}\left({\mathcal{A}}_{\mathrm{n}},,T_{\mathrm{n}}^{+},-,{\mathcal{A}}_{\mathrm{c}},,T_{\mathrm{c}}^{+}\right)\\ {\displaystyle }& {\displaystyle }+p_{\mathrm{n}\mathit{\alpha }}{\mathit{\phi }}_{\mathrm{n}-\mathrm{a}}\left({\mathcal{B}}_{\mathrm{n}},-,{\mathcal{B}}_{\mathrm{a}},+,{\mathcal{A}}_{\mathrm{n}},,T_{\mathrm{n}}^{+},-,{\mathcal{A}}_{\mathrm{a}},,T_{\mathrm{a}}^{+}\right)\\ {\displaystyle }& {\displaystyle }+\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right){\mathit{\phi }}_{\mathrm{n}-\mathrm{c}}{L}_{\mathrm{s}}\left[q_{\mathrm{sati}\mathrm{n}},,+,,{\displaystyle \frac{\partial q_{\mathrm{sati}\mathrm{n}}}{\partial T_{\mathrm{n}}}},\left(T_{\mathrm{n}}^{+},-,T_{\mathrm{c}}^{+}\right),,-,,q_{\mathrm{c}}^{+}\right]\\ {\displaystyle }& {\displaystyle }+p_{\mathrm{n}\mathit{\alpha }}{\mathit{\phi }}_{\mathrm{n}-\mathrm{a}}{L}_{\mathrm{s}}\left[q_{\mathrm{sati}\mathrm{n}},,+,,{\displaystyle \frac{\partial q_{\mathrm{sati}\mathrm{n}}}{\partial T_{\mathrm{n}}}},\left(T_{\mathrm{n}}^{+},-,T_{\mathrm{a}}^{+}\right),,-,,q_{\mathrm{a}}^{+}\right]\\ & {\displaystyle }& {\displaystyle }-{\mathrm{\Lambda }}_{\mathrm{g},\mathrm{1}}\left(T_{\mathrm{n},\mathrm{1}}^{+},-,T_{\mathrm{n},\mathrm{2}}^{+}\right)\mathit{\}}{p}_{\mathrm{ng}}.\end{array}$$ Note that Eq. () has been multiplied by $p_{\mathrm{ng}}$ since the snowpack must be made patch relative when solving the coupled equations. The $q_{\mathrm{sat}x}^{+}$ and longwave radiation terms have been linearized with respect to $T_x$ (the longwave radiation derivatives are given by Eq. ). The superscript $+$ corresponds to the values of variables at time $t+\mathrm{\Delta }t$, while the absence of a superscript indicates variables evaluated at time $t$. Note that we have defined ${\mathit{\phi }}_x={\mathit{\rho }}_{\mathrm{a}}/R_{\mathrm{a}x}$ (kg m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$) for simplicity. The thermodynamic variable, ${\mathcal{T}}_x$, in the sensible heat flux terms have been expressed as a function of $T_x$ using Eq. (). Several of the ${\mathcal{B}}_x$ terms have canceled out in the sensible heat flux terms in Eqs. ()–() since they are defined such that ${\mathcal{B}}_{\mathrm{c}}={\mathcal{B}}_{\mathrm{v}}={\mathcal{B}}_{\mathrm{g}}={\mathcal{B}}_{\mathrm{n}}$. Note that compared to Eqs. ()–(), the phase change terms (${\mathrm{\Phi }}_x$) do not appear in Eqs. ()–(). This is because they are evaluated as an adjustment after the energy budget and the fluxes have been computed.

In Eq. (), $T_{\mathrm{n},N_{\mathrm{n}}}^{\ast }$ represents a test temperature for the lowest snowpack layer. It is first computed using an implicit calculation of the combined snow–soil layers to get a first estimate of the snow–ground heat conduction inter-facial flux when simultaneously solving the surface energy budgets. The final snow temperature in this layer, $T_{\mathrm{n},N_{\mathrm{n}}}^{+}$, is computed afterwards within the snow scheme; any difference between the resulting conduction flux and the test flux in Eq. () is added to the soil as a correction at the end of the time step in order to conserve energy. In practice, this correction is generally small, especially since the snow fraction goes to unity very rapidly (i.e., for a fairly thin snowpack when using MEB; see Eq. ). Thus, in this general case, the difference between the test flux and the final flux arise only owing to updates to snow properties within the snow scheme during the time step. Since $T_{\mathrm{n},N_{\mathrm{n}}}^{\ast }$ is computed using an implicit solution method for the entire soil–snow continuum, it is also quite numerically stable. The use of a test flux permits a modular coupling between the snow scheme and the soil–vegetation parts of ISBA-MEB.

In order to solve Eqs. ()–() for the three unknown surface energy budget temperatures, $T_{\mathrm{v}}^{+}$, $T_{\mathrm{g},\mathrm{1}}^{+}$, and $T_{\mathrm{n},\mathrm{1}}^{+}$, equations for the six additional unknown surface energy budget temperatures, $T_{\mathrm{a}}^{+}$, $T_{\mathrm{c}}^{+}$, $q_{\mathrm{a}}^{+}$, $q_{\mathrm{c}}^{+}$, $T_{\mathrm{g},\mathrm{2}}^{+}$, and $T_{\mathrm{n},\mathrm{2}}^{+}$, must be defined. They can be expressed as linear equations in terms of $T_{\mathrm{v}}^{+}$, $T_{\mathrm{g},\mathrm{1}}^{+}$, and $T_{\mathrm{n},\mathrm{1}}^{+}$, and their derivations are presented in the remaining sections of this Appendix.

Atmospheric temperature and specific humidity

The first step in solving the surface energy budget is to eliminate the lowest atmospheric energy and water vapor variables from the snow surface energy budget equation. They will also be used to diagnose the final flux exchanges between the canopy air space and overlying atmosphere.

From Eq. (), the thermodynamic variable of the lowest atmospheric model variable at time $t+\mathrm{\Delta }t$ is defined as $${\mathcal{T}}_{N_{\mathrm{a}}}^{+}=B_{\mathcal{T},N_{\mathrm{a}}}+A_{\mathcal{T},N_{\mathrm{a}}}{H}^{+}.$$ Note that using Eq. (), we can rewrite Eq. () in terms of air temperature as $${T_{\mathrm{a}}}^{+}=B_{T_{\mathrm{a}}}+A_{T_{\mathrm{a}}}{H}^{+},$$ where $B_{T_{\mathrm{a}}}=\left(B_{\mathcal{T},N_{\mathrm{a}}},-,{\mathcal{B}}_{\mathrm{a}}\right)/{\mathcal{A}}_{\mathrm{a}}$, $A_{T_{\mathrm{a}}}=A_{\mathcal{T},N_{\mathrm{a}}}/{\mathcal{A}}_{\mathrm{a}}$, and $T_{\mathrm{a}}$ is shorthand for $T(k=N_{\mathrm{a}})$. Substitution of Eq. () for $H$ in Eq. () and solving for $T_{\mathrm{a}}^{+}$ yields $$T_{\mathrm{a}}^{+}={\stackrel{\mathrm{̃}}{B}}_{T_{\mathrm{a}}}+{\stackrel{\mathrm{̃}}{A}}_{T_{\mathrm{a}}}{T}_{\mathrm{c}}^{+}+{\stackrel{\mathrm{̃}}{C}}_{T_{\mathrm{a}}}{T}_{\mathrm{n}}^{+},$$ where $$\begin{array}{cc}& {\displaystyle }& {\displaystyle }C={\mathcal{A}}_{\mathrm{a}}\mathit{\{}\mathrm{1}+A_{T_{\mathrm{a}}}\left[{\mathit{\phi }}_{\mathrm{c}-\mathrm{a}},\left(\mathrm{1},-,p_{\mathrm{ng}},,p_{\mathit{\alpha }\mathrm{n}}\right),+,p_{\mathrm{ng}},,p_{\mathit{\alpha }\mathrm{n}},,{\mathit{\phi }}_{\mathrm{n}-\mathrm{a}}\right]\mathit{\}},& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{A}}_{T_{\mathrm{a}}}=A_{T_{\mathrm{a}}}{\mathit{\phi }}_{\mathrm{c}-\mathrm{a}}{\mathcal{A}}_{\mathrm{c}}\left(\mathrm{1},-,p_{\mathrm{ng}},,p_{\mathit{\alpha }\mathrm{n}}\right)/C,{\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{B}}_{T_{\mathrm{a}}}=\mathit{\{}{B}_{T_{\mathrm{a}}}-{\mathcal{B}}_{\mathrm{a}}+A_{T_{\mathrm{a}}}[\left(\mathrm{1},-,p_{\mathrm{ng}},,p_{\mathit{\alpha }\mathrm{n}}\right){\mathit{\phi }}_{\mathrm{c}-\mathrm{a}}\left({\mathcal{B}}_{\mathrm{c}},-,{\mathcal{B}}_{\mathrm{a}}\right)+\\ & {\displaystyle }& {\displaystyle }{p}_{\mathrm{ng}}{p}_{\mathit{\alpha }\mathrm{n}}{\mathit{\phi }}_{\mathrm{n}-\mathrm{a}}\left({\mathcal{B}}_{\mathrm{c}},-,{\mathcal{B}}_{\mathrm{a}}\right)\left]\mathit{\right\}}/C,& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{C}}_{T_{\mathrm{a}}}=A_{T_{\mathrm{a}}}{p}_{\mathrm{ng}}{p}_{\mathit{\alpha }\mathrm{n}}{\mathit{\phi }}_{\mathrm{n}-\mathrm{a}}{\mathcal{A}}_{\mathrm{c}}/C.\end{array}$$

In analogous fashion to determining the air temperature, the specific humidity of the lowest atmospheric model variable at time $t+\mathrm{\Delta }t$ is defined from Eq. () as $$q_{\mathrm{a}}^{+}=B_{q,\mathrm{a}}+A_{q,\mathrm{a}}{E}^{+},$$ where again the subscript $q$, a represents the values of the coefficients $A$ and $B$ for the lowest atmospheric model layer ($k=N_{\mathrm{a}}$). Substitution of Eq. () for $E$ in Eq. () and solving for $T_{\mathrm{a}}^{+}$ yields $$q_{\mathrm{a}}^{+}={\stackrel{\mathrm{̃}}{B}}_{q,\mathrm{a}}+{\stackrel{\mathrm{̃}}{A}}_{q,\mathrm{a}}{q}_{\mathrm{c}}^{+}+{\stackrel{\mathrm{̃}}{C}}_{q,\mathrm{a}}{q}_{\mathrm{sati}\mathrm{n}}^{+},$$ where the coefficients are defined as $$\begin{array}{}& {\displaystyle }& {\displaystyle }C=\mathrm{1}+A_{q,\mathrm{a}}\left[\left(\mathrm{1},,-,,p_{\mathrm{ng}},,p_{\mathit{\alpha }\mathrm{n}}\right),{\mathit{\phi }}_{\mathrm{c}-\mathrm{a}},,+,,{\mathit{\phi }}_{\mathrm{n}-\mathrm{a}},,h_{\mathrm{sn}},,p_{\mathit{\alpha }\mathrm{n}},,p_{\mathrm{ng}}\right],& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{A}}_{q,\mathrm{a}}=A_{q,\mathrm{a}}{\mathit{\phi }}_{\mathrm{c}-\mathrm{a}}\left(\mathrm{1},,-,,p_{\mathrm{ng}},,p_{\mathit{\alpha }\mathrm{n}}\right)/C,& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{B}}_{q,\mathrm{a}}=B_{q,\mathrm{a}}/C,& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{C}}_{q,\mathrm{a}}=A_{q,\mathrm{a}}{\mathit{\phi }}_{\mathrm{n}-\mathrm{a}}{h}_{\mathrm{sn}}{p}_{\mathit{\alpha }\mathrm{n}}{p}_{\mathrm{ng}}/C.\end{array}$$

Canopy air temperature and specific humidity

In order to close the energy budgets, $T_{\mathrm{c}}^{+}$ and $q_{\mathrm{c}}^{+}$ must be determined.

Assuming conservation of the heat flux between the different surfaces and the canopy air space, we have $$\begin{array}{cc}{\displaystyle }\left(\mathrm{1},-,p_{\mathrm{ng}},p_{\mathrm{n}\mathit{\alpha }}\right)H_{\mathrm{c}}^{+}& {\displaystyle }=p_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right)H_{\mathrm{n}-\mathrm{c}}^{+}\\ & {\displaystyle }& {\displaystyle }+\left(\mathrm{1},-,p_{\mathrm{ng}}\right)H_{\mathrm{g}}^{+}+H_{\mathrm{v}}^{+},\end{array}$$ which can be expanded as $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{\mathit{\phi }}_{\mathrm{c}-\mathrm{a}}\left(\mathrm{1},-,p_{\mathrm{ng}},,p_{\mathit{\alpha }\mathrm{n}}\right)\times \\ {\displaystyle }& {\displaystyle }\left({\mathcal{B}}_{\mathrm{c}},+,{\mathcal{A}}_{\mathrm{c}},,T_{\mathrm{c}}^{+},-,{\mathcal{B}}_{\mathrm{a}},-,{\mathcal{A}}_{\mathrm{a}},,T_{\mathrm{a}}^{+}\right)=\\ {\displaystyle }& {\displaystyle }{\mathcal{A}}_{\mathrm{c}}[{\mathit{\phi }}_{\mathrm{g}}\left(T_{\mathrm{g}}^{+},,-,,T_{\mathrm{c}}^{+}\right)\left(\mathrm{1},-,p_{\mathrm{ng}}\right)\\ {\displaystyle }& {\displaystyle }+{\mathit{\phi }}_{\mathrm{v}}\left(T_{\mathrm{v}}^{+},,-,,T_{\mathrm{c}}^{+}\right)\\ & {\displaystyle }& {\displaystyle }{\mathit{\phi }}_{\mathrm{n}-\mathrm{c}}\left(T_{\mathrm{n}}^{+},,-,,T_{\mathrm{c}}^{+}\right)p_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathit{\alpha }\mathrm{n}}\right)].\end{array}$$

Note that the above conservation equation does not include the part of the snow sensible heat flux, which is in direct contact with the atmosphere ($H_{\mathrm{n}-\mathrm{a}}$), since it was already accounted for in the expression for $T_{\mathrm{a}}^{+}$ via Eq. (). Eliminating $T_{\mathrm{a}}^{+}$ using Eq.  and solving for $T_{\mathrm{c}}^{+}$ yields $$T_{\mathrm{c}}^{+}=a_{T_{\mathrm{c}}}+b_{T_{\mathrm{c}}}{T}_{\mathrm{v}}^{+}+c_{T_{\mathrm{c}}}{T}_{\mathrm{g}}^{+}+d_{T_{\mathrm{c}}}{T}_{\mathrm{n}}^{+}$$ with the coefficients $$\begin{array}{cc}{\displaystyle }C=& {\displaystyle }{\mathit{\phi }}_{\mathrm{c}-\mathrm{a}}\left(\mathrm{1},,-,,p_{\mathrm{ng}},,p_{\mathit{\alpha }\mathrm{n}}\right)\left({\mathcal{A}}_{\mathrm{c}},-,{\mathcal{A}}_{\mathrm{a}},,{\stackrel{\mathrm{̃}}{A}}_{T\mathrm{a}}\right)+\\ & {\displaystyle }& {\displaystyle }{\mathcal{A}}_{\mathrm{c}}\left[{\mathit{\phi }}_{\mathrm{v}},+,{\mathit{\phi }}_{\mathrm{g}},,\left(\mathrm{1},-,p_{\mathrm{ng}}\right),+,{\mathit{\phi }}_{\mathrm{n}-\mathrm{c}},,p_{\mathrm{ng}},\left(\mathrm{1},-,p_{\mathit{\alpha }\mathrm{n}}\right)\right],& {\displaystyle }{a}_{T_{\mathrm{c}}}=& {\displaystyle }\left[{\mathit{\phi }}_{\mathrm{c}-\mathrm{a}}\left(\mathrm{1},,-,,p_{\mathrm{ng}},,p_{\mathit{\alpha }\mathrm{n}}\right)\left({\mathcal{B}}_{\mathrm{a}},-,{\mathcal{B}}_{\mathrm{c}},+,{\mathcal{A}}_{\mathrm{a}},,{\stackrel{\mathrm{̃}}{B}}_{T\mathrm{a}}\right)\right]/C,& {\displaystyle }{b}_{T_{\mathrm{c}}}=& {\displaystyle }{\mathcal{A}}_{\mathrm{c}}{\mathit{\phi }}_{\mathrm{v}}/C,& {\displaystyle }{c}_{T_{\mathrm{c}}}=& {\displaystyle }{\mathcal{A}}_{\mathrm{c}}{\mathit{\phi }}_{\mathrm{g}}\left(\mathrm{1},-,p_{\mathrm{ng}}\right)/C,{\displaystyle }{d}_{T_{\mathrm{c}}}=& {\displaystyle }\left[{\mathcal{A}}_{\mathrm{c}},,{\mathit{\phi }}_{\mathrm{n}-\mathrm{c}},,p_{\mathrm{ng}},\left(\mathrm{1},-,p_{\mathit{\alpha }\mathrm{n}}\right)\right.\\ & {\displaystyle }& {\displaystyle }\left..,+,,{\mathcal{A}}_{\mathrm{a}},,{\stackrel{\mathrm{̃}}{C}}_{T\mathrm{a}},,{\mathit{\phi }}_{\mathrm{c}-\mathrm{a}},\left(\mathrm{1},,-,,p_{\mathrm{ng}},,p_{\mathit{\alpha }\mathrm{n}}\right)\right]/C.\end{array}$$

In an analogous fashion for canopy air temperature determination, assuming conservation of the vapor flux between the different surfaces and the canopy air space, $$\begin{array}{cc}{\displaystyle }\left(\mathrm{1},-,p_{\mathrm{ng}},p_{\mathrm{n}\mathit{\alpha }}\right)E_{\mathrm{c}}^{+}& {\displaystyle }=p_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right)E_{\mathrm{n}-\mathrm{c}}^{+}\\ & {\displaystyle }& {\displaystyle }+\left(\mathrm{1},-,p_{\mathrm{ng}}\right)E_{\mathrm{g}}^{+}+E_{\mathrm{v}}^{+},\end{array}$$ which can be expanded using the definitions of the evaporative fluxes $E_x$ from Eqs. ()–() together with the definitions of $q_{\mathrm{g}}$ from Eq. () and $q_{\mathrm{a}}^{+}$ from Eq. () as $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{\mathit{\phi }}_{\mathrm{c}-\mathrm{a}}\left(\mathrm{1},-,p_{\mathrm{ng}},,p_{\mathit{\alpha }\mathrm{n}}\right)\times \\ {\displaystyle }& {\displaystyle }\left[q_{\mathrm{c}}^{+},\left(\mathrm{1},-,{\stackrel{\mathrm{̃}}{A}}_{q,\mathrm{a}}\right),,-,,{\stackrel{\mathrm{̃}}{B}}_{q,\mathrm{a}},,-,,{\stackrel{\mathrm{̃}}{C}}_{q,\mathrm{a}},,q_{\mathrm{sati}\mathrm{n}}^{+}\right]=\\ {\displaystyle }& {\displaystyle }[{\mathit{\phi }}_{\mathrm{g}}\left(h_{\mathrm{sg}},,q_{\mathrm{sat}\mathrm{g}}^{+},,-,,h_{\mathrm{a}},,q_{\mathrm{c}}^{+}\right)\left(\mathrm{1},-,p_{\mathrm{ng}}\right)+{\mathit{\phi }}_{\mathrm{v}}{h}_{\mathrm{sv}}\left(q_{\mathrm{sat}\mathrm{v}}^{+},,-,,q_{\mathrm{c}}^{+}\right)\\ & {\displaystyle }& {\displaystyle }{\mathit{\phi }}_{\mathrm{n}-\mathrm{c}}{h}_{\mathrm{sn}}\left(q_{\mathrm{sati}\mathrm{n}}^{+},,-,,q_{\mathrm{c}}^{+}\right)p_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathit{\alpha }\mathrm{n}}\right)].\end{array}$$ Owing to the linearization of the $q_{\mathrm{sat}x}$, terms about $T_x$, Eq. () can be solved for $q_{\mathrm{c}}^{+}$ as a function of the surface energy budget temperatures as $$q_{\mathrm{c}}^{+}=a_{q_{\mathrm{c}}}+b_{q_{\mathrm{c}}}{T}_{\mathrm{v}}^{+}+c_{q_{\mathrm{c}}}{T}_{\mathrm{g}}^{+}+d_{q_{\mathrm{c}}}{T}_{\mathrm{n}}^{+},$$ where the coefficients are defined as $$\begin{array}{cc}{\displaystyle }C=& {\displaystyle }{\mathit{\phi }}_{\mathrm{c}-\mathrm{a}}\left(\mathrm{1},,-,,p_{\mathrm{ng}},,p_{\mathrm{n}\mathit{\alpha }}\right)\left(\mathrm{1},-,{\stackrel{\mathrm{̃}}{A}}_{q,\mathrm{a}}\right)+{\mathit{\phi }}_{\mathrm{g}}{h}_N\left(\mathrm{1},-,p_{\mathrm{ng}}\right)\\ & {\displaystyle }& {\displaystyle }+{\mathit{\phi }}_{\mathrm{v}}{h}_{\mathrm{sv}}+{\mathit{\phi }}_{\mathrm{n}-\mathrm{c}}{h}_{\mathrm{sn}}{p}_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right),{\displaystyle }{a}_{q_{\mathrm{c}}}=& {\displaystyle }\mathit{\{}\left(\mathrm{1},,-,,p_{\mathrm{ng}},,p_{\mathrm{n}\mathit{\alpha }}\right){\mathit{\phi }}_{\mathrm{c}-\mathrm{a}}{\stackrel{\mathrm{̃}}{B}}_{q,\mathrm{a}}+{\mathit{\phi }}_{\mathrm{v}}{h}_{\mathrm{sv}}\left(q_{\mathrm{sat}\mathrm{v}},,-,,{\displaystyle \frac{\partial q_{\mathrm{sat}\mathrm{v}}}{\partial T_{\mathrm{v}}}},,T_{\mathrm{v}}\right)\\ {\displaystyle }& {\displaystyle }+{\mathit{\phi }}_{\mathrm{g}}{h}_{\mathrm{sg}}\left(q_{\mathrm{sat}\mathrm{g}},,-,,{\displaystyle \frac{\partial q_{\mathrm{sat}\mathrm{g}}}{\partial T_{\mathrm{g}}}},,T_{\mathrm{g}}\right)\left(\mathrm{1},-,p_{\mathrm{ng}}\right)\\ & {\displaystyle }& {\displaystyle }+{\mathit{\phi }}_{\mathrm{n}-\mathrm{c}}{h}_{\mathrm{sn}}\left(q_{\mathrm{sati}\mathrm{n}},,-,,{\displaystyle \frac{\partial q_{\mathrm{sati}\mathrm{n}}}{\partial T_{\mathrm{n}}}},,T_{\mathrm{n}}\right)& {\displaystyle }& {\displaystyle }{p}_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right)\mathit{\}}/C,& {\displaystyle }{b}_{q_{\mathrm{c}}}=& {\displaystyle }{h}_{\mathrm{sv}}{\mathit{\phi }}_{\mathrm{v}}{\displaystyle \frac{\partial q_{\mathrm{sat}\mathrm{v}}}{\partial T_{\mathrm{v}}}}/C,& {\displaystyle }{c}_{q_{\mathrm{c}}}=& {\displaystyle }{h}_{\mathrm{sg}}{\mathit{\phi }}_{\mathrm{g}}{\displaystyle \frac{\partial q_{\mathrm{sat}\mathrm{g}}}{\partial T_{\mathrm{g}}}}\left(\mathrm{1},-,p_{\mathrm{ng}}\right)/C,& {\displaystyle }{d}_{q_{\mathrm{c}}}=& {\displaystyle }{h}_{\mathrm{sn}}{\mathit{\phi }}_{\mathrm{n}-\mathrm{c}}{\displaystyle \frac{\partial q_{\mathrm{sati}\mathrm{n}}}{\partial T_{\mathrm{n}}}}{p}_{\mathrm{ng}}\left(\mathrm{1},-,p_{\mathrm{n}\mathit{\alpha }}\right)/C.\end{array}$$

Sub-surface temperatures

The sub-surface conduction heat fluxes (Eqs. –) can be expressed in compact form as $$G_{x,k}^{+}={\mathrm{\Lambda }}_{x,k}\left(T_{x,k}^{+},-,T_{x,k+\mathrm{1}}^{+}\right),$$ where ${\mathrm{\Lambda }}_{x,k}$ represents the ratio of the inter-facial thermal conductivity to the thickness between the mid-points of contiguous layers ($k$ and $k+\mathrm{1}$). Using the methodology described in Appendix  for the atmospheric diffusion scheme, the soil and snow heat diffusion equation (both using the form of Eq. ) can be defined in an analogous fashion as $$T_{\mathrm{g},k}^{+}=B_{\mathrm{g},k}+A_{\mathrm{g},k}{T}_{\mathrm{g},k-\mathrm{1}}^{+}\left(k,=,\mathrm{2},,,N_{\mathrm{g}}\right),$$ where the coefficients $B_{\mathrm{g},k}$ and $A_{\mathrm{g},k}$ are determined during the upward sweep (first step of the tri-diagonal solution) from the base of the soil to the sub-surface soil and snow layers as described by . The resulting coefficients for the soil are defined as $$\begin{array}{cc}& {\displaystyle }C=& {\displaystyle }\left({\mathcal{C}}_{\mathrm{g}k},/,\mathrm{\Delta },t\right)+{\mathrm{\Lambda }}_{\mathrm{g}k-\mathrm{1}}+{\mathrm{\Lambda }}_{\mathrm{g}k}\left(\mathrm{1},-,A_{\mathrm{g}k+\mathrm{1}}\right),{\displaystyle }{B}_{\mathrm{g}\mathrm{i}}=& {\displaystyle }\left[\left({\mathcal{C}}_{\mathrm{g}k},/,\mathrm{\Delta },t\right),,T,\mathrm{g},,k,,+,,{\mathrm{\Lambda }}_{\mathrm{g}k},,B_{\mathrm{g}k+\mathrm{1}}\right]\\ & {\displaystyle }& {\displaystyle }/C\left(\mathrm{2},\le ,k,\le ,N_{\mathrm{g}},-,\mathrm{1}\right),& {\displaystyle }{A}_{\mathrm{g}k}=& {\displaystyle }{\mathrm{\Lambda }}_{\mathrm{g}k-\mathrm{1}}/C.\end{array}$$ The same form holds for the snow layers. The upward sweep is performed before the evaluation of the energy budget; thus, Eq. () is used to eliminate $T_{\mathrm{g},\mathrm{2}}^{+}$ and $T_{\mathrm{n},\mathrm{2}}^{+}$ from Eqs. () and (), respectively. To do this, the sub-surface implicit fluxes in Eqs. () and () can be expressed, respectively, as $$\begin{array}{}& {\displaystyle }{G}_{\mathrm{g},\mathrm{1}}^{+}=& {\displaystyle }{\mathrm{\Lambda }}_{\mathrm{g},\mathrm{1}}\left[T_{\mathrm{g},\mathrm{1}}^{+},\left(\mathrm{1},,-,,A_{\mathrm{g},\mathrm{2}}\right),,+,,B_{\mathrm{g},\mathrm{2}}\right],& {\displaystyle }{G}_{\mathrm{n},\mathrm{1}}^{+}=& {\displaystyle }{\mathrm{\Lambda }}_{\mathrm{n},\mathrm{1}}\left[T_{\mathrm{n},\mathrm{1}}^{+},\left(\mathrm{1},,-,,A_{\mathrm{n},\mathrm{2}}\right),,+,,B_{\mathrm{n},\mathrm{2}}\right].\end{array}$$

Surface stresses

Using the same surface–atmosphere coupling methodology as for temperature and specific humidity, the $u$ wind component in the lowest atmospheric model layer can be expressed as $$u_{\mathrm{a}}^{+}=B_{u\mathrm{a}}+A_{u\mathrm{a}}{\mathit{\tau }}_x^{+}.$$ The surface $u$ component momentum exchange with the atmosphere is expressed as $${\mathit{\tau }}_x^{+}=-u_{\mathrm{a}}^{+}\left[\left(\mathrm{1},-,p_{\mathrm{ng}},,p_{\mathrm{n}\mathit{\alpha }}\right),{\mathit{\phi }}_{D\mathrm{c}-\mathrm{a}},,+,,p_{\mathrm{ng}},,p_{\mathrm{n}\mathit{\alpha }},{\mathit{\phi }}_{D\mathrm{n}-\mathrm{a}}\right],$$ where it includes stresses from the snow-buried and non-snow-buried portions of the surface consistent with the fluxes of heat and water vapor. For simplicity, we have defined $${\mathit{\phi }}_{{\mathrm{D}}_x}={\mathit{\rho }}_{\mathrm{a}}{V}_{\mathrm{a}}{C}_{{\mathrm{D}}_x}$$ and $C_{\mathrm{D}}$ is the surface drag coefficient, which is defined following . Eliminating ${\mathit{\tau }}_x^{+}$ from Eq. () using Eq. () gives $$u_{\mathrm{a}}^{+}={\displaystyle \frac{B_{u\mathrm{a}}}{\mathrm{1}+A_{u\mathrm{a}}{\mathit{\phi }}_{D_{\mathrm{c}}}}},$$ where for convenience we have defined the average drag coefficient as $${\mathit{\phi }}_{D_{\mathrm{c}}}=\left(\mathrm{1},-,p_{\mathrm{ng}},,p_{\mathrm{n}\mathit{\alpha }}\right){\mathit{\phi }}_{D\mathrm{c}-\mathrm{a}}+p_{\mathrm{ng}}{p}_{\mathrm{n}\mathit{\alpha }}{\mathit{\phi }}_{D\mathrm{n}-\mathrm{a}}.$$ The net $u$-momentum flux from the surface to the canopy air space is expressed as $${\mathit{\tau }}_x^{+}=-{\displaystyle \frac{B_{u\mathrm{a}}{\mathit{\phi }}_{D_{\mathrm{c}}}}{\left(\mathrm{1},,+,,A_{u\mathrm{a}},,{\mathit{\phi }}_{D_{\mathrm{c}}}\right)}}.$$ Finally, the scalar friction velocity can be computed from $$u^{\ast }={\left({\displaystyle \frac{{\mathit{\phi }}_{D_{\mathrm{c}}}{V}_{\mathrm{a}}^{+}}{{\mathit{\rho }}_{\mathrm{a}}}}\right)}^{\mathrm{1}/\mathrm{2}},$$ where $V_{\mathrm{a}}^{+}$ is the updated wind speed (computed from $u_{\mathrm{a}}^{+}$ and $v_{\mathrm{a}}^{+}$). Note that $v_{\mathrm{a}}^{+}$ and ${\mathit{\tau }}_y^{+}$ are computed in the same manner, but using $B_{\mathrm{v}\mathrm{a}}$ from the atmosphere (note that $A_{\mathrm{v}\mathrm{a}}=A_{u\mathrm{a}}$).

Summary: final solution of the implicitly coupled equations

The fully implicit solution of the surface and atmospheric variables proceeds for each model time step as follows:

• Within the atmospheric model, perform the downward sweep of the tri-diagonal matrix within the turbulent diffusion scheme of the atmospheric model to obtain the $A_{\mathit{\varphi },k}$ and $B_{\mathit{\varphi },k}$ coefficients for each diffused variable ($\mathit{\varphi }=\mathcal{T}$, $q$, $u$, and $v$) for each layer of the atmosphere ($k=\mathrm{1},N_{\mathrm{a}}$). Update ${\mathcal{A}}_{\mathrm{a}}$ and ${\mathcal{B}}_{\mathrm{a}}$, then pass these values along with the aforementioned coupling coefficients at the lowest atmospheric model layer (i.e., $A_{T,\mathrm{a}}$, $B_{T,\mathrm{a}}$, $A_{q,\mathrm{a}}$, $B_{q,\mathrm{a}}$, $A_{u,\mathrm{a}}$, $B_{u,\mathrm{a}}$, and $B_{\mathrm{v},\mathrm{a}}$) to the land surface model. These coefficients are then used to eliminate $T_{\mathrm{a}}^{+}$ and $q_{\mathrm{a}}^{+}$ from the implicit surface energy budget equations (Eqs. –).

• Within the land surface model, perform the upward sweep of the tri-diagonal matrix within the soil and snow layers to determine the $A_{\mathrm{n},k}$, $B_{\mathrm{n},k}$, $A_{\mathrm{g},k}$, and $B_{\mathrm{g},k}$, coefficients for the soil and snow layers (from soil-layer $N_{\mathrm{g}}$ to layer 2, and again from soil-layer $N_{\mathrm{g}}$ to layer 2 of the snow scheme). Note that coefficients for layer 1 of the snow and soil schemes are not needed since they correspond to the linearized surface energy budgets (next step).

• Within the land surface model, the expressions for $T_{\mathrm{a}}^{+}$ (Eq. ), $q_{\mathrm{a}}^{+}$ (Eq. ), $T_{\mathrm{c}}^{+}$ (Eq. ), $q_{\mathrm{c}}^{+}$ (Eq. ), $T_{\mathrm{g},\mathrm{2}}^{+}$ (Eq. )and $T_{\mathrm{n},\mathrm{2}}^{+}$ (Eq. ) can now be substituted into the energy budget equations (Eqs. –), which can then be readily solved for $T_{\mathrm{v}}^{+}$, $T_{\mathrm{g},\mathrm{1}}^{+}$, and $T_{\mathrm{n},\mathrm{1}}^{+}$.

• Within the land surface model, perform back substitution (using $T_{\mathrm{g},\mathrm{1}}^{+}$ as the upper boundary condition) to obtain $T_{\mathrm{g},k}^{+}$ for soil layers $k=\mathrm{2},N_{\mathrm{g}}$ using Eq. ().

• Within the land surface model, call the explicit snow-process scheme to update the snow scheme temperature, $T_{\mathrm{n},k}^{+}$, and the snow mass variables for snow layers $k=\mathrm{2},N_{\mathrm{n}}$. The implicit snow surface fluxes, $R_{\mathrm{n},\mathrm{n}}^{+}$, $H_{\mathrm{n}}^{+}$ and $E_{\mathrm{n}}^{+}$, are used as the upper boundary condition along with the implicit soil temperature, $T_{\mathrm{g},\mathrm{1}}^{+}$, to compute the updated lower snowpack boundary condition (i.e., the snow–soil inter-facial flux, $G_{\mathrm{gn}}$).

• Within the land surface model, compute $V_{\mathrm{a}}^{+}$ (see Sect. ). Diagnose $T_{\mathrm{a}}^{+}$, $T_{\mathrm{c}}+$, $q_{\mathrm{a}}^{+}$ and $q_{\mathrm{c}}^{+}$ (again, using the equations mentioned in step 3) in order to compute the updated (implicit) fluxes. The updated evapotranspiration (Eqs. –) and snowmelt water mass fluxes are used within the hydrology schemes to update the different water storage variables for the soil and vegetation canopy (Eqs. –).

• Within the atmospheric model, perform back substitution (using $H^{+}$, $E^{+}$, ${\mathit{\tau }}_x^{+}$ and ${\mathit{\tau }}_y^{+}$ as the lower boundary conditions: Eq. ) to obtain updated profiles (or turbulent tendencies, depending on the setup of the atmospheric model) of ${\mathcal{T}}_k$, $q_k$, $u_k$ and $v_k$ for atmospheric layers $k=\mathrm{1},N_{\mathrm{a}}$. Finally, the updated upwelling shortwave, $\mathrm{SW}\uparrow$, and implicit longwave flux, $\mathrm{LW}{\uparrow }^{+}$ (or equivalently, the effective emissivity and implicit longwave radiative temperature, $T_{\mathrm{rad}}^{+}$) are returned to the atmospheric model as lower boundary conditions for the respective radiative schemes.

Alternately, in offline mode, $A_{\mathit{\varphi },\mathrm{a}}=\mathrm{0}$ and $B_{\mathit{\varphi },\mathrm{a}}={\mathit{\varphi }}_{\mathrm{a}}$ in the driving routine in step 1, and the solution procedure ends at step 6. Finally, if multiple patches and/or tiles are being used within the grid call of interest, the corresponding fractional-area-weighted fluxes are passed to the atmospheric model in step 7.

Acknowledgements

This work was initiated within the international HIRLAM consortium as part of the ongoing collective effort to improve the SURFEX platform for research and operational hydrological and meteorological applications. We wish to make a posthumous acknowledgement of the contribution to this work by Joël Noilhan, who was one of the original supporters of this project and helped initiate this endeavor; his scientific vision was essential at the early stages of this work. We wish to thank other contributors to this development in terms of discussions and evaluation, such as G. Boulet, E. Martin, J.-C. Calvet, P. Le Moigne, C. Canac, and G. Aouad. The technical support of S. Faroux of the SURFEX team is also greatly appreciated. We also wish to thank E. Lebas for preparing the MEB schematic. Part of this work was supported by a grant from Météo-France.Edited by: S. Valcke Reviewed by: two anonymous referees