ISBA-MEB (SURFEX v8.1): model snow evaluation for local-scale forest sites

Accurate modeling of the effect of snow cover on the surface energy and mass fluxes is required from land surface models. The Interactions between Soil–Biosphere–Atmosphere (ISBA) model uses a composite soil–vegetation approach that has limitations when representing snow and soil phase change processes in areas of high vegetation cover since it does not explicitly represent the snowpack lying on the ground below the canopy. In particular, previous studies using ISBA have pointed out that the snowpack ablation tends to occur to early in the season in forest regions in the Northern Hemisphere. The multi-energy balance (MEB) version of ISBA has been developed recently, to a large degree, to address this issue. A vegetation layer, which is distinct from the soil, has been added to ISBA and new processes are now explicitly represented, such as snow interception and an understory litter layer. To evaluate the behavior of this new scheme in a cold forested region, long-term offline simulations have been performed for the three BERMS forest sites located in Saskatchewan, Canada. It is shown that the new scheme leads to an improved energy budget representation, especially in terms of the ground and sensible heat fluxes, with decreases in root-mean-square error (RMSE) of 77 % and 18 %, respectively. A positive impact for soil temperatures, consistent with the improvement of the ground heat flux, is obtained, particularly in terms of bias, which is reduced from −6.2 to −0.1 K at a 10 cm soil depth on average for the three sites and 12 studied years. The impact of using MEB on the snowpack simulation is a better agreement with observations during the snow season, especially concerning the last day of snow in the season: errors are on the order of 1 d averaged over the three sites and all of the years using MEB, which represents a reduction in error of 20 d compared to the composite scheme. The analysis shows that this improvement is mostly caused by the ability of MEB to represent a snowpack that nearly completely covers the soil below the canopy and that decouples the soil from the atmosphere, while keeping a close coupling between the vegetation and the atmosphere.

1 Introduction

Forests cover approximately one-third of world's land surface area, one-third of which consists of boreal forests in subarctic and cold continental climates. In these regions, snowpack can last more than half of the year and can modify the surface roughness and thermal and radiative properties, and thereby it has a significant impact on the fluxes of momentum, heat, and water mass between the surface and the atmosphere or the soil. Vegetation canopy processes in forests modulate the behavior (accumulation and melting) of the snowpack on the ground. Notably, snowfall can be intercepted by the canopy leaves and branches, where it can be sublimated or melted before unloading to the forest floor (Pomeroy et al., 1998; Storck et al., 2002; Bartlett et al., 2006). In addition, downwelling shortwave radiative fluxes are attenuated by the sheltering effect of the canopy (Harding and Pomeroy, 1996), while the longwave radiation reaching the below-canopy snow surface is generally enhanced compared to its atmospheric component due to longwave radiation emission by the canopy and trunks (Gouttevin et al., 2015; Todt et al., 2018). The snowpack constitutes a very efficient thermal insulating material that decreases the cooling of the soil compared to a snow-free surface (Zhang, 2005; Grundstein et al., 2005), which in turn can have a significant impact on soil freezing and thawing and thus on the permafrost depth (Stieglitz et al., 2003; Paquin and Sushama, 2015).

Land surface models (LSMs) seek to provide realistic simulations of snow evolution, which implies that they have the ability to represent the previously mentioned first-order processes. An accurate representation of the impact of snow cover on the energy and water balances is required for land surface reanalysis products, particularly in cold regions (Carrera et al., 2015) and for operational regional-scale hydrological modeling for which snowmelt is a key driver of discharge (e.g. Habets et al., 2008; Snow et al., 2016). In addition, more physically based multi-layer snow schemes have been developed for operational numerical weather prediction (Dutra et al., 2010), including explicit forest canopy formulations (Yang et al., 2011), and such schemes have also been developed for climate modeling (e.g. Oleson et al., 2010; Decharme et al., 2016).

The consensus of global climate model (GCM) predictions for the current century is that high-latitude regions will continue to warm at an accelerated rate compared to other regions of the globe, in large part owing to the positive snow albedo feedback (Flanner et al., 2011; Qu and Hall, 2014). This mechanism is considered to be a driver of the observed Arctic amplification of the current global warming (e.g. Bony et al., 2006; Chapin et al., 2005; Serreze and Barry, 2011). But it is known that the spread in surface albedo feedback among different CMIP5 GCMs is particularly large in the boreal forest zone (Qu and Hall, 2014), and this is in large part owing to the representation of snow masking by vegetation (Thackeray et al., 2015). Thackeray et al. (2018) state that the main reason for the spread in surface albedo is owing to structural aspects of the LSMs (the representation of snowpack, the vegetation canopy and their interactions) rather than the parameter values used by these schemes.

In the 1990s, a series of Model Inter-comparison Projects (MIPs) were initiated in order to intercompare and evaluate the LSM state-of-the-art representation of cold season processes with the goal of determining which aspects of the schemes were affecting performance and causing model spread and also to provide guidance for future model developments. Multiple MIPs at the local scale, for which detailed measurements of snow processes exist, have been done over the past 20 years, such as the Programme for Intercomparison of Land-Surface Parameterization Schemes (PILPS) Phase 2d (Slater et al., 2001) and SnowMIP Phase 1 (Etchevers et al., 2004) and Phase 2 (Essery et al., 2009; Rutter et al., 2009). The PILPS-Phase 2e experiment (Bowling et al., 2003) looked at the combined effect of multiple cold-season processes at the regional scale over a Scandinavian catchment. The Rhone-AGGregation MIP (Boone et al., 2004) evaluated the snow depth simulations of an ensemble of LSMs at numerous observation sites in the French Alps: they determined that liquid water retention (owing to a storage capacity in the snow and also possible refreezing of this liquid water) was a key processes required for simulating the accurate timing and amount of snowmelt and thus discharge in a high alpine catchment. All of the aforementioned MIPs used observation-based forcing as boundary conditions to the LSMs in offline (decoupled from the atmosphere) mode. More recently, the Earth System Model Snow Model Intercomparison Project (ESM-SnowMIP: Krinner et al., 2018) extended the intercomparison to the global scale and also used fully coupled GCM-LSM models. Note that in particular, SnowMIP2 and ESM-SnowMIP evaluations highlighted the difficulties LSMs have in modeling snow in forested sites compared to open sites for a large number of LSMs.

In order to represent certain snow processes in forested regions (e.g., the shielding effect of the canopy, the downwelling longwave enhancement, or the fractional coverage of snow lying on the ground), many LSMs have adopted an explicit representation of the vegetation canopy as opposed to the composite schemes that consider the ground and the canopy as a unique entity. In the late 1980s and early 1990s, the so-called two-source energy budget method began to be implemented into GCMs. In this approach, the surface (which can consist in soil, snow, or a blend thereof) is distinct from the overlying bulk vegetation canopy, each computing their own fluxes and having explicit parameters. The first and most simplified version was proposed by Deardorff (1978). Based on that approach, Sellers et al. (1986) proposed one of the first comprehensive schemes for use in a GCM that inspired and still resembles many of the two-source LSMs in use today. Over time, more variations of this type of approach have emerged, such as LSMs using simplified treatments for certain processes, such as radiative transfer and a significantly reduced number of input parameters (more easily adapted for use in GCMs), while still retaining the overall explicit canopy (Xue et al., 1991). Some LSMs have further split the canopy into two layers (Saux-Picart et al., 2009) representing the overstory and understory vegetation layers, or by separating the trunk from the vegetation canopy with a focus on including important longwave radiative impacts that can be critical for below-canopy snowpack evolution (Gouttevin et al., 2012). More recently, models have been developed using a multi-layer vegetation canopy for GCM applications (Ryder et al., 2016) with improvements to the explicit treatment of turbulence (Bonan et al., 2018).

The Interaction Soil–Biosphere–Atmosphère (ISBA) LSM is part of the platform SURFEX (SURFace EXternalisée: Externalized Surface) software platform being developed at Météo-France in collaboration with multiple international partners (Masson et al., 2013). SURFEX is used in operational systems such as numerical weather prediction within the global atmospheric model ARPEGE (Action de Recherche Petite Echelle Grande Echelle) operational at Météo France or the limited-area model AROME (Seity et al., 2011) or hydrological and land surface analysis systems (Habets et al., 2008). It is also used within the CNRM-CM6 climate model (Decharme et al., 2019), which is participating in the Coupled Model Intercomparison Project phase 6 (CMIP6) project (Eyring et al., 2016). Several key improvements have recently been made to ISBA that impact the representation of cold season processes: soil water phase changes are governed using the Gibbs free energy concept, liquid water and temperature computations now extend over more soil layers (resulting in a higher vertical resolution in the upper layers, along with the ability to represent soil temperature for very deep soil layers), the explicit snowpack includes more layers enabling high resolution at both the atmosphere and soil interfaces and the thermal conductivity and albedo parameterizations have been improved, and finally the soil can include soil organics (Decharme et al., 2016). Despite these improvements, ISBA still has difficulties simulating the snowpack evolution and soil temperatures for forested sites as evidenced by SnowMIP2. Along the same lines, Decharme et al. (2019) identified a precocious snowmelt over boreal forest regions by a global simulation which led to a springtime peak of river discharge over all Arctic basins that was too early. This issue is related to the conceptual aspect of the composite model: the snowpack cannot be represented between the upper ground layer and the canopy, and thus a compromise has to be found, notably for the snow fraction calculation, which results in a partial snow coverage in forests at all times, thus enabling an excessive direct coupling between the atmosphere and both the below-canopy ground and snowpack. The main effect is to cool the snowpack and especially the ground too much below dense forests.

For these reasons, the Multi Energy Balance (MEB, Boone et al., 2017) option was recently implemented within ISBA to allow an explicit and distinct representation of the upper ground and vegetation layers. Radiative transfer models for shortwave and longwave radiation are improved, interception of snow by the canopy layer is added and turbulent fluxes formulation are adapted to the new design. A parameterization has also been added to model the litter on the ground (Napoly et al., 2017), which reduces soil evaporation and heat exchanges with the soil. The MEB option has been evaluated on a large number of local forest sites (Napoly et al., 2017), but little attention has been paid specifically to its impact on the simulation of snow until now.

In the current study, the ability of ISBA-MEB to model the snowpack is evaluated using data from the Boreal Ecosystem Research Study (BERMS), which covers a 12-year (1 January 1999–31 December 2010) observational period for three distinct (aspen-, jack pine- and black spruce-dominated) Canadian forest sites (Bartlett et al., 2006). The motivation behind the selection of these sites is to study areas for which MEB has the dramatic impact in terms of reducing model bias. The current operational ISBA (single composite soil–vegetation scheme) is used as a reference model, and thus the aim of this study is to highlight the performance of the new MEB option for the modeling of the snowpack and the related variables. After a presentation of the main contrasting characteristics of each model and a description of the forest sites, the two options are evaluated and compared using the data from the BERMS sites with a focus on the modeling of certain key features of the snowpack (e.g., the timing and length of the melting period, the snow energy budget, and the impact on soil temperatures). Several sensitivity tests are also summarized, focusing on the most uncertain parameters of the new MEB option.

2 Model

The ISBA model is developed within the SURFEX platform (SURface EXternalisée, Masson et al., 2013), and version 8.1 is used in this study. There are multiple parameterization options available, notably those which govern soil thermal and hydrological fluxes, snowpack physics, and the new explicit vegetation canopy and forest litter options. Note that ISBA within SURFEX includes the notion of explicit sub-grid patches, which represent different types of plant functional types or land classes explicitly, but in the current study we use a single-patch representation (thus, throughout the text, we refer to patch or grid cell interchangeably). The options most pertinent to the current study are described below.

2.1 ISBA: default configuration

ISBA uses a so-called composite approach which is defined herein as using a single energy balance for the combined soil–vegetation surface Fig. 1a. The properties of soil and vegetation are aggregated depending on the fraction that vegetation occupies (veg) in the considered grid cell (Noilhan and Planton, 1989; Noilhan and Mahfouf, 1996). The snow fraction is also used in the aggregation of properties of the surface. Its formulation, which is detailed in Sect. 2.1.4 and plotted on Fig. 2, results from a compromise between an atmospheric and soil point of view and represent one of the significant limitation of the model. Since the inception of ISBA, many developments have been made to improve the representation of physical processes as the knowledge of key processes, the quality and spatial and temporal coverage and resolution of input datasets, and the computing speed have improved. In this paper, we use the default ISBA configuration which is defined for the current study as follows.

• The soil water and energy transfers are simulated using the diffusive approach option (DIF) (Boone et al., 2000; Decharme et al., 2011) that uses multiple (here 12) layers to solve the Fourier and Darcy laws throughout the soil. The soil parameters are derived from soil texture using pedotransfer functions based on Clapp and Hornberger (1978) classification. The impact of soil organic carbon (SOC) on thermal and hydraulic properties is also used (Decharme et al., 2016).

• The parameterization of the stomatal resistance used to calculate the forest transpiration models the functional coupling between the stomatal resistance and the net assimilation of CO₂ (Ag-s, Calvet et al., 1998). An option to simulate the evolution of the leaf area index (LAI) prognostically is not activated in the current study since estimated values are available and thus imposed.

• The snowpack is modeled using a multi-layer physically based explicit snow option (ES) that partitions the layers based on snow depth (Boone and Etchevers, 2001). Since that time, multiple improvements over the ensuing 15 years have been implemented and are described by Decharme et al. (2016). The key physical processes are briefly summarized in Sect. 2.1.3.

Figure 1 Schematics of ISBA (a) and MEB (b) during a snow period on a forest site.

In the following section, we describe the aspects related to snow representation in the model that differ between the default version of ISBA and the new MEB option.

2.1.1 Energy budget

The energy budget equations for the composite surface soil–vegetation (hereafter simply referred to as the composite layer) and upper snow layer are expressed as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\mathcal{C}}_{\mathrm{s}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{s}}}{\mathrm{d}t}}=& \left(R_{\text{net\hspace{0.17em}s}}-{\mathrm{LE}}_{\mathrm{s}}\right.\left.-H_{\mathrm{s}}\right)\left(\mathrm{1}-p_{\mathrm{sn}}\right)+G_{\mathrm{ns}}{p}_{\mathrm{sn}}\\ & -G_{\mathrm{s},\mathrm{1}}+L_{\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{s},\mathrm{1}},\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\mathcal{C}}_{\mathrm{n},\mathrm{1}}{\displaystyle \frac{\partial T_{\mathrm{n},\mathrm{1}}}{\partial t}}=& R_{\text{net\hspace{0.17em}n}}-H_{\mathrm{n}}-{\mathrm{LE}}_{\mathrm{n}}-{\mathit{\tau }}_{\mathrm{n},\mathrm{1}}{\mathrm{SW}}_{\text{net\hspace{0.17em}n}}\\ & +{\mathit{\xi }}_{\mathrm{n},\mathrm{1}}-G_{\mathrm{n},\mathrm{1}}+L_{\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{n},\mathrm{1}},\end{array}\end{array}$$

where T_s (K) is the temperature of the composite surface and T_n,1 (K) represents the temperature of the uppermost layer of snow. 𝒞_s and 𝒞_n,1 (J K⁻¹ m⁻²) are the effective heat capacities of the composite and upper snow layers, respectively. Both budgets use a relatively thin layer (for soil and snow: on the order of several centimeters at the maximum) in order to be able to properly model the surface temperature diurnal cycle. R_net s and R_net n (W m⁻²) correspond to the net radiative fluxes for the soil and the snowpack, respectively. In the same way, LE_s, LE_n, H_s, and H_n (W m⁻²) are the latent and sensible heat fluxes. G_s,1 and G_n,1 (W m⁻²) are the conductive fluxes from the composite and snow surface layers to the corresponding sub-surface layers, respectively. The conductive heat flux between the base of the snowpack and the composite layer is represented by G_ns. The effective heating (or cooling) rate 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) is represented by ξ_n,1. Note that the integral of ξ_n over the entire snowpack depth is zero at the end of each time step. The phase change terms (defined following the convention of freezing minus melting, expressed in kg m⁻² s⁻¹) are represented by Φ_s and Φ_n, respectively, and L_f represents the latent heat of fusion (J kg⁻¹).

The fraction of the soil–vegetation surface covered by snow is p_sn, thus the surface soil layer is in contact simultaneously with both the base of the snowpack and the atmosphere when p_sn<1 (which represents a critical difference with MEB which will be discussed further in a subsequent section). Also note that the budget in Eq. (2) is snow-relative: in order to obtain the total energy budget and the net fluxes for a patch containing snow, all of the terms in the snow energy budget are multiplied by p_sn and then added to Eq. (1). Several of the terms most critical to cold season processes in Eqs. (1)–(2) are described in more detail in the following sections.

2.1.2 Radiative transfer

The surface net radiation of the soil–vegetation and snowpack are given by

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{R}_{\mathrm{net}}=& p_{\mathrm{sn}}\left({\mathrm{SW}}_{\text{net\hspace{0.17em}n}}+{\mathrm{LW}}_{\text{net\hspace{0.17em}n}}\right),\\ & +\left(\mathrm{1}-p_{\mathrm{sn}}\right)\left({\mathrm{SW}}_{\text{net\hspace{0.17em}s}}+{\mathrm{LW}}_{\text{net\hspace{0.17em}s}}\right),\end{array}\end{array}$$

where SW and LW represent the shortwave and longwave radiative flux components, respectively. Part of the incoming shortwave radiation received by the snowpack is transmitted through the uppermost snow layer, and this energy loss is expressed as τ_n,1 SW_net n, where τ is a dimensionless transmission coefficient, where the snow surface net shortwave radiation is

$$\begin{array}{}\text{(4)}& {\mathrm{SW}}_{\text{net\hspace{0.17em}n}}={\mathit{\alpha }}_{\mathrm{n}}\mathrm{SW}\downarrow ,\end{array}$$

where SW↓ is the atmospheric downwelling shortwave radiation. The transmission function is described in detail in Decharme et al. (2016). The total surface net shortwave radiation is defined using a so-called composite albedo defined as

$$\begin{array}{}\text{(5)}& {\mathit{\alpha }}_{\mathrm{s}}=\mathrm{veg}{\mathit{\alpha }}_{\mathrm{v}}+(\mathrm{1}-\mathrm{veg}){\mathit{\alpha }}_{\mathrm{g}},\end{array}$$

The total surface effective albedo of the snow soil-vegetation composite surface (α_eff) is then defined by weighting the contribution of each surface:

$$\begin{array}{}\text{(6)}& {\mathit{\alpha }}_{\mathrm{eff}}=p_{\mathrm{sn}}{\mathit{\alpha }}_{\mathrm{n}}+\left(\mathrm{1}-p_{\mathrm{sn}}\right){\mathit{\alpha }}_{\mathrm{s}},\end{array}$$

with α_n, α_v, and α_g the snow, vegetation, and ground albedos, respectively. Note that no explicit shortwave transmission through the canopy is modeled.

The net longwave radiation for either surface is defined as

$$\begin{array}{}\text{(7)}& {\mathrm{LW}}_{\mathrm{net}X}={\mathit{ϵ}}_X\left(\mathrm{LW}\downarrow -\mathit{\sigma }{T}_{X,\mathrm{1}}^{\mathrm{4}}\right),\end{array}$$

where X represents either s or n, σ is the Stefan–Boltzmann constant, LW↓ is the downwelling atmospheric radiation, and ϵ represents the emissivity. The effective emissivity (ϵ_eff) of the surface is then defined in a fashion analogous to the effective albedo as

$$\begin{array}{}\text{(8)}& {\mathit{ϵ}}_{\mathrm{eff}}=p_{\mathrm{sn}}{\mathit{ϵ}}_{\mathrm{n}}+(\mathrm{1}-p_{\mathrm{sn}})[\mathrm{veg}{\mathit{ϵ}}_{\mathrm{v}}+(\mathrm{1}-\mathrm{veg}\left){\mathit{ϵ}}_{\mathrm{g}}\right],\end{array}$$

with ϵ_n, ϵ_v and ϵ_g representing the snow, vegetation, and ground emissivity, respectively. This effective emissivity is then used to compute the effective surface radiative temperature (from the explicitly computed upwelling longwave fluxes from the snow and composite surfaces), which is required by the longwave radiative scheme when coupled to an atmospheric model.

2.1.3 Snow processes

The snowpack model (ISBA-ES) is a multi-layer snow model of intermediate complexity (Boone and Etchevers, 2001; Decharme et al., 2016). The model current uses a default of 12 layers to model the physical processes involved in the snowpack such as solar energy absorption, compaction, snowmelt, water percolation, and refreezing of meltwater. The snow albedo is based on a snow historical variable and considers up to three spectral bands. Readers are referred to the aforementioned references for more details.

2.1.4 The snow fraction

In the ISBA composite method, the effective fraction of the grid cell covered by snow (p_sn) is the average between the fraction of snow covering the vegetation and the one covering the ground. It is calculated as follows:

$$\begin{array}{}\text{(9)}& {\displaystyle }{p}_{\mathrm{sn}}& {\displaystyle }=\mathrm{veg}{p}_{\mathrm{snv}}+(\mathrm{1}-\mathrm{veg})p_{\mathrm{sng}},\text{(10)}& {\displaystyle }{p}_{\mathrm{snv}}& {\displaystyle }=\min \left(\mathrm{1.0},{\displaystyle \frac{D}{D+\mathrm{2}{z}_{\mathrm{0}}}}\right),\text{(11)}& {\displaystyle }{p}_{\mathrm{sng}}& {\displaystyle }=\min \left(\mathrm{1.0},{\displaystyle \frac{D}{D_{\mathrm{g}}}}\right),\end{array}$$

where the p_snv and p_sng values correspond to snow fraction over the vegetation and the ground, respectively, and D is the total snow depth (m). Note that several options for the parameterizations of p_snv and p_sng exist in SURFEX; however, for the current study, Eqs. (10)–(11) represent the default used with the multi-layer soil and snow schemes. There is no explicit canopy snow reservoir, thus only a masking effect of the vegetation cover is modeled. In order to avoid excessive bare-soil evaporation, the default value of the veg parameter is 0.95 for forests (used herein), z₀ (m) corresponds to the surface roughness, which is calculated as 0.13 times the vegetation height and D_g (m) is a snow depth threshold set to 0.01 m (the default value). As a result, for a forest patch, the maximum value of p_sn reaches a maximum of approximately 0.2 for a forest height of 11 m (corresponding to one of the sites in the current study) as shown in Fig. 2. This implies that part of the soil–vegetation composite surface is always in direct contact with the atmosphere regardless of the snow depth, thus the insulating effect of the snowpack is reduced in a forest compared to a bare or low-vegetation covered surface using the composite option. The reasoning for a parameterization resulting in a low p_snv value over forests represents a compromise between insulating the soil surface while not burying a forest, which would result in an unrealistic coupling with the overlying atmosphere, notably in terms of the total upwelling shortwave radiative flux.

Figure 2 Snow fraction for different values of SWE for ISBA and MEB over a forest of 11 m height and assuming a density of 200 kg m⁻³.

2.2 ISBA-MEB: explicit vegetation canopy

The ISBA-MEB option treats up to three fully coupled distinct surface energy budgets (Fig. 1b) which are the snow surface, the bulk vegetation canopy, and the ground, which is characterized in the current study as a litter layer (Napoly et al., 2017). The reader is referred to Boone et al. (2017) for an extended description of the various assumptions of the MEB approach, its full set of governing equations, and its numerical aspects. Compared to the classic ISBA approach, there are two additional prognostic heat storage variables, which are the vegetation temperature, T_v, and the litter temperature, T_L. There are also three new hydrological prognostic variables: the snow liquid water equivalent intercepted by the vegetation canopy, W_rn (kg m⁻²) and the liquid and liquid water equivalent ice stored in the litter layer, W_l (kg m⁻²) and W_li (kg m⁻²), respectively. In addition, note that the vegetation fraction parameter used to aggregate soil and vegetation properties in the composite method, veg, is not used in ISBA-MEB since the canopy and soil properties are modeled explicitly. Running the model with the MEB option adds an extra cost of 6 %.

2.2.1 Energy budget

The MEB coupled energy budget equations includes an additional energy budget for the bulk vegetation canopy, and in the current study, the additional litter energy budget equation is also included, which results in a modified upper boundary condition for the uppermost soil temperature. The new and modified energy budget equations are as follows:

$$\begin{array}{}\text{(12)}& {\displaystyle }{\mathcal{C}}_{\mathrm{v}}{\displaystyle \frac{\partial T_{\mathrm{v}}}{\partial t}}=R_{\text{net\hspace{0.17em}v}}-H_{\mathrm{v}}-{\mathrm{LE}}_{\mathrm{v}}+L_{\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{v}},\text{(13)}& {\displaystyle }\begin{array}{rl}{\mathcal{C}}_{\mathrm{l}}{\displaystyle \frac{\partial T_{\mathrm{l}}}{\partial t}}& =(\mathrm{1}-p_{\mathrm{sng}})(R_{\text{net\hspace{0.17em}l}}-H_{\mathrm{l}}-{\mathrm{LE}}_{\mathrm{l}})\\ & +p_{\mathrm{sng}}(G_{\text{nl}}+{\mathit{\tau }}_{\mathrm{n},N_{\mathrm{n}}}{\mathrm{SW}}_{\text{net\hspace{0.17em}n}})-G_{\mathrm{l}}+L_{\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{l}},\end{array}\end{array}$$

where T_v and T_l are the temperatures (K) of the bulk vegetation and litter layers, respectively, while 𝒞_v and 𝒞_l correspond to the effective heat capacities (J K⁻¹ m⁻²). R_net v, R_net l, H_v, H_l, LE_v and LE_l (W m⁻²) represent the same quantities as in Eqs. (1)–(2) but for the bulk vegetation and litter layers. Note that 𝒞_v includes the heat capacities of intercepted solid and liquid water. Note that the snow surface energy budget equation (Eq. 2) is unchanged, however, the definition of the net radiation term has changed. G_nl (W m⁻²) is the conductive heat flux between the lowest snow and the litter layer, and G_l (W m⁻²) is the conductive flux between the litter and the uppermost ground layer. Thus, in contrast to ISBA, the uppermost soil temperature in MEB is only modulated by conductive heat flux divergence and phase changes as follows:

$$\begin{array}{}\text{(14)}& {\mathcal{C}}_{\mathrm{g},\mathrm{1}}{\displaystyle \frac{\partial T_{\mathrm{g},\mathrm{1}}}{\partial t}}=G_{\mathrm{l}}-G_{\mathrm{g},\mathrm{1}}+L_{\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{g},\mathrm{1}},\end{array}$$

where 𝒞_g,1 represents the surface soil heat capacity (i.e., with no vegetation effects included). As in the equation of Sect. 2.1.1, water phase change terms, Φ_v and Φ_l, are included for the vegetation and the litter respectively.

2.2.2 Radiative transfer through the canopy

MEB represents the explicit radiative transfer through the vegetation for shortwave and longwave fluxes using classical approaches, and it is fully described in Sect. 2.4.2 of Boone et al. (2017). A few key aspects that are pertinent to the current study are described herein. The model uses the classic representation of the canopy as plane parallel surface with a canopy absorption defined as follows:

$$\begin{array}{}\text{(15)}& {\mathit{\sigma }}_{\mathrm{LW}}=\mathrm{1}-\exp (-{\mathit{\tau }}_{\mathrm{LW}}\mathrm{LAI}),\end{array}$$

where LAI corresponds to the leaf area index (m² m⁻²) and τ_LW is a coefficient which is set to 0.4 as a default. The model results can be impacted by this parameter, and some sensitivity tests are presented in Sect. 4.2.1. Note that compared to the composite scheme, MEB increases the downwelling longwave radiation (towards the soil and snowpack) and thus the below-canopy net longwave radiation by including an emission from the canopy (which can be significantly larger than the atmospheric component in cold or dry climates, such as in the current study).

The shortwave radiative transfer scheme is described in Carrer et al. (2013). It uses an explicit multi-layer computation that accounts for different characteristics of the vegetation, such as the leaf area index, the clumping index, direct and diffuse radiation components, the thickness of the leaves, and the zenith angle. The main outputs are the bulk-canopy reflected, transmitted, and absorbed radiation components, and the corresponding photosynthetically active radiation (PAR) used within the photosynthesis scheme. Compared to the composite scheme, the main impacts are that the downwelling radiation at the surface ground and snowpack is attenuated mainly as a function of the LAI. In addition, because the snowpack is generally below the canopy (for the forest heights considered in the current study), the upwelling shortwave radiation is generally significantly reduced compared to the composite scheme since the forest can effectively mask the surface. In the current study, the reflected shortwave radiation is merely a diagnostic, but in a coupled atmospheric model, the shortwave exchange can be significantly modified (the total upwelling shortwave radiation can be significantly reduced in Boreal forest zones, but exploring this impact is beyond the scope of the current study).

2.2.3 New and modified snow processes

Both the composite and MEB options are coupled to the identical version of the ES snow scheme. The impact of using the MEB option on snow processes compared to the composite version can be briefly summarized as follows.

• Only the snow fraction over the ground is considered (veg=0 in MEB) so that the snow fraction is simply defined as p_sn=p_sng. This implies that p_sn is generally much closer to unity for MEB than for ISBA in forests (based on Eqs. 9–11 and the discussion in Sect. 2.1.4). This implies a greater coverage of the ground by snow in MEB, while the canopy is totally exposed to the overlying atmosphere in MEB. Note that when snow depth becomes comparable with the height of the vegetation (for example, for shrubs or grasses), another parameter described in Boone et al. (2017) is introduced. However, it is not relevant in the current study for the forest heights considered herein.

• An explicit canopy snow reservoir is considered in MEB, which includes interception, unloading, and both freezing of intercepted liquid water and melting of snow. It is based on Hedstrom and Pomeroy (1998) and the implementation in MEB is described in detail in Boone et al. (2017).

• It is currently assumed that the impact of intercepted snow on the total canopy albedo is negligible. This is based on the results of Pomeroy and Dion (1996). They indicated that the scattering and multiple reflections of light due to intercepted snow, combined with the high probability for the reflected light to reach the underside of an overlying branch and leaves, implies that trees actually behave as light traps. They concluded that intercepted snow has no significant impact on the canopy shortwave albedo or on the net radiative exchange. Some models, such as ECHAMS, (Roesch et al., 2001), use a simple parametrization to consider the impact of intercepted snow on the total effective albedo. It mainly consists of considering an effective canopy albedo affected by the amount of snow in snow leaf reservoir.

• The fluxes from the snowpack are calculated using the specific humidity and temperature of the so called “canopy” air space (Fig. 1b) instead of the forcing “air” layer when using ISBA (Fig. 1a). This permits some feedback between the surface and the atmosphere on the fluxes.

• The below-canopy wind speed that impacts the ground-based snowpack is reduced owing to the attenuation of the wind speed due to vegetation that would tend to reduce sublimation, however, the fractional coverage is generally considerably larger, and thus generally snowpack sublimation is increased. In addition, sublimation can also occur from intercepted snow.

3 Data

The BERMS program (Boreal Ecosystem Research and Monitoring Sites) sites are used in this study. The studied period ranges from 1 January 1999 at 00:00 to 31 December 2010 at 23:30 UTC, corresponding to 12 years of measurements. The three sites are located in Saskatchewan, Canada, and are described in detail in Bartlett et al. (2006). Their distinguishing characteristics are listed in Table 1 and can be briefly summarized as follows.

• OAS. This site is dominated by 21 m average-height old aspen that naturally regenerated after a fire in 1919. A 2 m high understory composed mainly of hazelnut is present. The ground is characterized by (from the surface downward) an 8–10 cm layer of forest litter, a peat layer, and finally a sandy clay loam soil.

• OBS. Old black spruce is the dominant tree species of this site. Trees have an average height of 12 m. The understory is comprised of shrubs and herbs, mosses and lichens that are situated on sandy loam and sandy soil.

• OJP. The old jack pine site is approximately 14 m high and is composed of a very sparse understory (alder, bearberry, cranberry, and lichens) over a coarse sandy soil.

Table 1 Main physical characteristics of the sites.

The full set of meteorological observations needed to force an LSM (downwelling all wavelength solar and atmospheric radiation fluxes, air temperature and humidity, pressure, liquid and solid precipitation, and wind speed above the forest canopy) is available at half-hourly time steps over the full period, along with data which enable a detailed description of the vegetation characteristics (such as LAI, albedo; see Table 2).

Table 2 Model characteristics prescribed for the three different sites.

In order to evaluate the model performance, measurements of turbulent fluxes, upwelling short and longwave radiation, soil temperature, and volumetric water content profiles are also available at a 30 min time step. Snow depth was measured every 30 min during the duration of the ground-based snowpack. In addition, manual measurements of snow water equivalent were made up to six times per year. The observed shortwave radiation being transmitted through the canopy (reaching the ground or snowpack surface: SW_g) was derived from the photosynthetically active radiation (PAR). The data were filtered to account for measurement error due to a direct flux on the sensor around midday, which caused very high peaks. The filter is based on the surrounding three points, using the following threshold:

$$\begin{array}{}\text{(16)}& \mathit{\delta }=\mathrm{abs}\left[{\mathrm{SW}}_{\mathrm{g}}\right[i]-\mathrm{0.5}({\mathrm{PSW}}_{\mathrm{g}}[i-\mathrm{1}]+{\mathrm{SW}}_{\mathrm{g}}[i+\mathrm{1}]\left)\right]\end{array}$$

and thus we reject the observation at i if δ≥100 W m². Note that this threshold is somewhat arbitrary because the anomalous peaks are quite large relative to the surrounding values and generally last for one time step. Due to the lack of a frost and snow cleaning system on the PAR sensor, we did not use measurements corresponding to these conditions in the evaluation.

The energy balance closure for these sites has been calculated as follows:

$$\begin{array}{}\text{(17)}& \text{closure}={\displaystyle \frac{\stackrel{\mathrm{‾}}{H+\mathrm{LE}}}{\stackrel{\mathrm{‾}}{R_{\mathrm{net}}-G-S}}},\end{array}$$

where the overbar corresponds to averages over the study period. The storage and ground heat fluxes are not considered here as they were not measured. In addition, it was assumed that they are, on average, negligible compared to the net radiation when averaged over such a long period. The energy balance closure was 84 %, 91 %, and 90 % for the OBS, OJP, and OAS sites, respectively. This closure is deemed to be satisfactory for the analysis in the current study, especially with respect to the study of Wilson et al. (2002), which found an average closure of 80 % over the Fluxnet network sites.

4 Results

All of the simulations are one-dimensional, and we assume that the forests are fairly homogeneous around the observation sites. Options for the explicit multi-layer vertical soil heat and water transfer (DIF) and ground-based snowpack (ES) are used, along with the Ag-s stomatal resistance formulation (described in Sect. 2), and thus only the impact of the MEB option is evaluated. In the following, we will refer to the different experiments as follows: MEB for the experiment using the new Multi Energy Balance option and ISBA for the default experiment. To evaluate the new option, a statistical analysis is performed that is based on simulated fluxes, soil variables, and the snowpack characteristics. Then, the study focuses on some specific periods where snow plays a key role governing the surface and sub-surface processes. Finally, a sensitivity analysis is performed to test several new MEB parameters that are the most likely to influence the snow processes. Model input parameters have been chosen that correspond to site measurements where possible. For the remaining parameters, we use the physiographic database developed for SURFEX (ECOCLIMAP, Champeaux et al., 2005) and the soil parameters from the HWSD dataset (harmonized world soil database, Nachtergaele and Batjes, 2012). The main parameters are given in Table 2.

4.1 Evaluation

4.1.1 Energy fluxes

One of the most critical fluxes in coupled land–atmosphere simulations over cold regions is the upwelling shortwave radiation, SW↑. The simulated flux is relatively close between the two experiments and to the observations, as shown in Figs. 3b and Fig. 4d, e, f, with averaged root-mean-square errors (RMSEs) over all sites and years of 7.7 and 7.1 W m⁻² for MEB and ISBA, respectively (Table 3). Improving the modeling of the reflected solar radiation would mostly consist of improving the quality of the input parameters, i.e. albedos (visible and near-infrared values for the soil and the vegetation) and LAI. At the deciduous OAS site in winter, the LAI is low (about 1.0 m² m⁻²) and the SW↑ is overestimated, notably in MEB. We suspect that a stem area index (SAI) should be explicitly considered to lower the effect of snow below the canopy on the effective albedo, especially on such a forest of 22 m height, consistent with results from Napoly et al. (2017). The solar radiation that passes through the canopy is only modeled with the new MEB option. When data were available, the simulation of the radiation that is transmitted through the canopy is rather well modeled (Fig. 3a). Unfortunately, the quality of the data was not sufficient enough when LAI was low at the OAS deciduous site to confirm the assumption of the importance of including a SAI. In winter, solar radiation remains relatively low at this site and barely affects the surface energy balance, and thus this issue is not addressed for the moment.

Table 3 RMSE for the ISBA-MEB and ISBA experiments for fluxes SW↑, LW↑, H, LE, and G calculated over half-hourly data.

Figure 3 Composite of monthly diurnal cycle at the OJP site. MEB is in red, ISBA in blue, measurements are indicated by a solid line, and adjusted measurements are represented using a dashed black line. Adjusted measurements respect the energy balance closure and are calculated following the method of Twine et al. (2000). As a visual aid, the area between the latter two is shaded and model outputs should fall within this area.

The impact of MEB for the OBS site is the opposite to that at the OAS site. At this site, the LAI is relatively large (3.5 m² m⁻²; see Table 1). Thus, in MEB, the total effective surface albedo is approximately equal to the canopy albedo. In ISBA, there is always a fraction of snow visible to the overlaying atmosphere. Even though it is relatively low (from 10 %–20 %, as shown in Fig. 2), it can result in an overestimation of the total reflected shortwave radiation, especially when the snow is fresh and the snow albedo is relatively large compared to that of the vegetation. This effect is seen in Fig. 4e. The ISBA bias arises mainly from an over-estimation of the effective surface albedo in February and March (not shown).

Sensible heat flux (H) is well simulated over the three sites with MEB compared to ISBA as shown in Figs. 3d and 4g, h, i, with an average RMSE of 48.4 W m⁻² (respectively 58.9 W m⁻²) and average bias of 4.1 W m⁻² (respectively −1.0 W m⁻²). This result is consistent with Napoly et al. (2017), who showed that the large overestimation of the ground heat flux diurnal amplitude from ISBA, confirmed in this study (Fig. 3f), results in a lack of energy in turbulent fluxes and mostly in H as LE is notably limited by the evaporative demand. This overestimation is largely decreased with MEB due to the shielding effect of the canopy (thereby reducing the net solar radiation at the below-canopy surface) and the insulating effect of the explicit litter layer reduces the heat exchanges with the below-surface soil layers. During periods with snow cover, this improvement is even more marked due to the presence of the snowpack and is investigated in the next section.

Figure 4 Scatterplots for the different fluxes and the three sites using only observations with snow on the ground. MEB is in red, and ISBA is in blue.

Simulations are also improved for the latent heat flux (LE) with an average RMSE of 37.1 W m⁻² for MEB and 47.3 W m⁻² for ISBA, and an average BIAS of 6.6 W m⁻² for MEB and 9.9 W m⁻² for ISBA. The main differences appear during spring when ISBA tends to overestimate the total evapotranspiration mainly owing to an excessive soil evaporation (despite the fact that only a 5 % soil fraction is prescribed (RMSE and bias are shown in Tables 3 and 4, respectively). The default veg value of 0.95 for forests has been tuned to avoid excess bare-soil evaporation in ISBA. In MEB, no tuned parameter is required to limit bare-soil evaporation and it is generally lower than in ISBA owing to a lower surface roughness length (since the surface roughness of the soil is fixed to a few centimeters at most in MEB but can be several tens of centimeters in ISBA since the soil and vegetation properties are aggregated) and diminished wind speeds owing to the frictional effects of the canopy. In addition, the total annual bare-soil evaporation is further reduced in MEB over the seasonal cycle since the ablation is later. Finally, the litter layer also has an impact on reducing the ground evaporation, and further explanation will be given in the next section.

Table 4 BIAS for the ISBA-MEB and ISBA experiments for fluxes SW↑, LW↑, H, LE, and G calculated over half-hourly data.

The sublimation of snow represents 27 % (12 % from the snowpack itself and 15 % from the intercepted snow by the canopy) of the total snowfall using MEB, whereas it is only 2 % with ISBA. This change occurs for essentially two reasons: (i) the snow fraction parameterization gives a low value of snow cover for ISBA compared to MEB (Fig. 2), which weights the fluxes, and (ii) with MEB the interception of the snow by the canopy is explicitly considered and allows more sublimation. Even if no observations can confirm these differences, studies have estimated that in forests sublimation might represent several tens of percent of the annual snowfall (Pomeroy and Dion, 1996) and may exceed 30 % (Montesi et al., 2004). More recently, Molotch et al. (2007) measured a ratio between sublimation of the snowpack and total sublimation of snow of 0.45 for a forest in Colorado (at 3000 m), which is quite close to the values of 0.42, 0.45, and 0.51 found here for the OJP, OBS, and OAS sites, respectively.

4.1.2 Snow

The 12-year average annual cycle of the snowpack evolution is shown in Fig. 5 for the ISBA (blue curve) and default MEB (solid red line) simulations. The statistical scores calculated on the last day of snow when comparing simulation to measurements are shown in Table 5. The last day of snow was defined using the two following conditions: (i) the first time when total snow depth (SND) <0.02 m was identified and (ii) where the average SND over the ensuing 2 weeks remained below this threshold value. These simple criteria were found to determine the timing of the melt of the snowpack quite accurately without being mistaken with a possible late snow event. In addition, this result is not very sensitive to the chosen threshold value: the idea is simply to eliminate short-term snow cover events occurring after the main ablation. The average and standard deviation of the BIAS between modeled and observed last day of snow are shown in Table 5. With the ISBA option, snow melts on average 24 d too early. Using MEB leads to an improvement in the simulation of the fluxes shown in the previous section as well as the snowpack depth, with an average RMSE of 7.7 cm and BIAS of 0.4 cm (Table 8). The most important effect appears in springtime when MEB simulates ablation later, with an average BIAS in last day of snow of only 1 d (too early) averaged over the 3 sites and full time period. In order to better illustrate the ability of the model to represent the snowpack below the canopy, Fig. 6 shows the SND and SWE evolution over fairly a representative period consisting of 3 consecutive years (early 2001 to early 2004). The relatively less frequently observed values of the snow water equivalent SWE allow a confirmation of the good representation of the timing of snow melt. In addition, Fig. 6 seems to indicate that the snow density is well modeled since underestimation or overestimation of SND and SWE are consistent for both models. However, SWE measurements are not numerous enough to accurately confirm this good behavior. Over the whole period of study, ISBA has an RMSE of 13.8 cm and BIAS of −4.0 cm (Table 8). Errors mainly arise because the melt of the snowpack occurs too early in the spring season, and this behavior is consistent for all years studied. Figure 7 displays different parameters at the OJP site from 25 to 31 March 2004, which correspond to a melting period. The significant overestimation of the ISBA surface soil temperature fluctuations is obvious (Fig. 7c) as it almost perfectly follows the temperature measured at 5 m above the soil. Thus, it leads a large conductive heat flux between the soil and the snowpack (Fig. 7d) (on the order of several hundred W m⁻²), which is unrealistic compared to net radiation (Fig. 7e). This is explained by the relatively low fraction (10 %; see Eq. 9) occupied by the snowpack for a forest in the composite model ISBA. This fraction allows the model to simulate a rather good effective total albedo (Fig. 3c), but as a result approximately 90 % of soil is not shielded by snow and is strongly coupled to the atmospheric forcing.

Table 5 Average and standard deviation of the BIAS between model and observations of the last day of snow expressed in number of days. The OAS site has only 6 years of snow observations compared to 9 years for both the OBS and OJP.

Figure 5 Composite of snow depth annual (July to June) cycles for the three sites: MEB is in red, ISBA is in blue, and observations are in black. The dashed red curve corresponds to a version of MEB with no snow interception by the canopy.

Figure 6 Snow depth (a, c, e) and snow water equivalent (b, d, f) for the three sites from 1 July 2001 to 1 July 2004. MEB is in red, ISBA is in blue, and observations are in black.

As spring begins, the atmospheric temperature gets closer to 0 ^∘C and solar radiation starts to increase. With ISBA, the ground temperature can easily rise to over 0 ^∘C, as the heat capacity of that layer is low and part of the ground surface is directly exposed to the atmosphere (again owing to a relatively low p_sn value compared to MEB). Once the ground temperature exceeds 0 ^∘C, the conductive flux between the snow and the ground (Fig. 7d) is negative, indicating that the ground is warming up the snowpack from below. The early melt of the snowpack in ISBA is thus due, in large part, to the energy received from the combined ground–vegetation layer (Fig. 7c). In MEB, the insulation of the soil from the snowpack is total as the horizontal coverage of the snow is more realistic. The flux coming from the ground is very close to 0 W m⁻² (Fig. 7d). Thus, the melt of the snowpack comes almost entirely from above (as the snow becomes thin, some solar energy can warm the ground below the snowpack, thereby melting the snow from below using MEB as well, but this effect tends to be quite small compared to melting induced by the surface flux of heat into the snowpack). The net radiation (Fig. 7e) received by the first layer of the snowpack is higher in MEB than in ISBA due to the longwave enhancement effect, and this causes the snowpack to melt at a speed more comparable to the measurements (Fig. 7a).

Figure 7 Multiple variables at the OJP site from 25 to 31 March 2004, which corresponds to a melting period of the snowpack. MEB is in red, ISBA is in blue, and observations are in black. T_G1 is the temperature of the first layer of the surface (i.e. the composite for ISBA and the first layer of the ground for MEB).

A period before the ablation of the snowpack is shown in Fig. 8. Certain fluxes from the ISBA energy budget (H and G) are quite different compared to the observations. Indeed, because the snowpack does not cover the full grid, the available energy is used to warm up or cool down the surface soil temperature, which provokes a strong amplitude of G instead of being released to the atmosphere through H. With MEB, two prognostic temperatures (Fig. 1) are used, T_l for the surface litter, which barely varies in time, and T_v for the exchanges with the atmosphere, which is related to an explicit heat capacity of the vegetation (lower than the composite heat capacity of ISBA). These two temperatures, which are totally uncorrelated due to the snowpack that occupies the full surface of the ground, lead to a much improved modeling of energy fluxes.

Figure 8 Multiple variables at the OJP site from 25 January to 1 February 2008. MEB is in red, ISBA is in blue, and observations are in black. For T, the red curve is T_C (MEB), the dotted red curve is T_G1 (MEB), and the blue curve is T_G1 (ISBA).

4.1.3 Soil temperature and water content

The overestimation of the ground heat flux amplitude by ISBA impacts not only energy exchanges with the atmosphere through H but also the soil temperatures. With the direct contact of about 90 % of the composite layer with the atmosphere, the soil temperature at a depth of 10 cm calculated from ISBA can drop to below −20 ^∘C in winter months (Fig. 9), whereas observed temperatures at this depth are only slightly negative: this feature is common for the three sites and during the entire period (results are shown here for 3 years for ease of visual inspection). This leads to a significant temperature BIAS averaged over the three sites and full time period at 10 cm depth of −2.9 K and an RMSE of 6.8 K (Tables 6 and 7). Owing to the insulating effect of the snowpack, MEB is much closer to observations with an average BIAS of 0.1 K and RMSE of 2.0 K.

Figure 9 Soil water content and temperature at 7.5 and 10 cm depth, respectively, for the three sites from 2 July 2001 to 2 July 2004. MEB is in red, ISBA is in blue, and observations are in black. On the WG (soil water) graphs, the dotted lines represent the liquid and solid water.

Table 6 RMSE for the ISBA-MEB and ISBA experiments for soil temperature at 5, 20, and 100 cm calculated over half-hourly data.

Table 7 BIAS for the ISBA-MEB and ISBA experiments for soil temperature at 5, 20, and 100 cm calculated over half-hourly data.

In ISBA, the increased exposure to cold atmospheric conditions leads to a cold bias that extends to at least 1 m depth throughout the season (the 12-year average seasonal cycle of soil temperature for the three sites is shown in Fig. 10), which is the maximum depth of the soil temperature observations. This leads to significantly more soil freezing with depth during early winter. In spring, even after the snowpack ablation, the frozen water component remains significant in the deep soil layers. This cold bias is mainly owing to the under-estimated impact of the insulating effect of snow (a low p_sn in ISBA) Finally, note that even if the MEB-simulated soil temperatures warm a bit more slowly than the observations (as evidenced by the small delay and slight tilt of the annual temperature wave compared to the observations indicating a bit more inertia in MEB versus the observations), MEB provides a much improved soil temperature simulation.

Figure 10 Average annual contours of soil temperature between surface and 100 cm deep at the OJP site.

The near-surface (7.5 cm depth) modeled soil liquid water content (Fig. 9) agrees reasonably well with the observations for both versions of the model except for the OJP site. The overestimated values at this site are likely to be due to the definition of the soil characteristics, which are defined based on soil texture information. A noticeable difference between the two models is that the water content curves are generally more flat, with the MEB option in months outside of summer than ISBA, which is in better agreement with the observations. ISBA indeed occasionally melts the entire snowpack erroneously, as shown in Fig. 6, leading to short periods of ice melting and unrealistic peaks of liquid water content. The impact of changes in soil freezing between MEB and ISBA on drainage and runoff are thus expected in regional or global studies, and this could further have an impact on the hydrological cycle (notably river flow) in such regions. This issue using ISBA was identified by Decharme et al. (2019), who found a precocious springtime peak of river discharge over all Arctic basins. Thus, it is anticipated that MEB should at least improve this bias in future large-scale hydrological studies using SURFEX in both coupled and offline modes.

4.2 Sensitivity tests

Several sensitivity tests were performed and the results are summarized here. The analysis focuses on three parameters and one process for which the values are considered to be uncertain and for which the snowpack is potentially sensitive. For each parameter, values were tested for each site over a range (either based on the literature or physical reasoning) and compared to the default value defined in Boone et al. (2017). Statistical scores were calculated only when snow was observed on the ground.

4.2.1 Canopy longwave radiation transmission

The τ_LW parameter is an absorption coefficient that is used to calculate the LW radiation transmitted through the canopy (Eq. 15), and it weights the canopy emission to the soil and the atmosphere. The original default value of τ_LW is 0.5 (Boone et al., 2017) and values from 0.1 to 1.0 using increments of 0.1 are tested in the current study. This range covers values based on a literature survey as discussed in the aforementioned reference, although 1.0 is quite large and is tested simply for numerical reasons. As τ_LW increases, the canopy transmission decreases and the canopy emission increases (increasing the longwave radiation received by the snowpack). Figure 11 shows the RMSE calculated for each value of this parameter over the 12-year period for each site for the identified four most impacted state variables and flux. The sensitivity to this parameter is relatively high, notably for low values. For each variable, errors are lowest for τ_LW values in the range 0.3–0.4. For further increases in this parameter, the RMSE stabilizes for LW↑ and H, and starts to increase for SND and T_G. This behavior is consistent for all three sites (as shown). Thus, the value of 0.4, quite close to the default one, has been selected to be the new default value (and it has been used for the results presented in previous sections).

Figure 11 RMSE calculated for different values of τ_LW for each site during the snow period and for the sensible heat flux (a), snow depth (b), soil temperature at 2 cm (c), and LW↑ (d).

4.2.2 Litter thickness

The litter thickness has been identified as a key parameter of MEB (Napoly et al., 2017). Indeed, it affects the thermal and hydrological fluxes and state variables in the model since its thickness models the litter surface energy budget and water storage capacity. Its value can be very specific for a particular site and can evolve in time, and, in addition, values are hard to determine at a large scales. Its variation is generally in the range from 0.01 to 0.10 m based on a literature survey shown in Napoly et al. (2017). In addition, this range has been selected for essentially two additional reasons: the first is that model tests have shown that results degrade for thicknesses below 0.01 m since the assumption of the existence of a continuous litter layer becomes physically dubious and if it is too thin, not to mention numerical issues can arise since the surface energy budget is computed within this layer. Second, when the layer exceeds approximately 0.10 m, the diurnal cycle is highly damped (to levels which are unrealistic): a multi-layer litter model would be preferable within the MEB model structure in this case. These two issues are discussed in more detail in Napoly et al. (2017): based on the aforementioned study, the MEB default value for litter thickness is constant in time and set to 0.03 m. However, in the current study, tests showed no significant sensitivity to this parameter during snow periods (<10 % variation of RMSE of the tested values compared to the default value) on most of the state variables and total surface–atmosphere fluxes. Only the soil temperatures were found to be significantly impacted, with optimized RMSE values of 1–2 K obtained with litter thickness values at or above 0.06 m, instead of 3–4 K with 0.01 m. Note that these differences are essentially due to the initial state of the ground when the snow season begins since the litter effect is active during the entire year.

4.2.3 Roughness length for heat and water vapor

The ratio of the vegetation roughness lengths of momentum to heat (and vapor), $r_{z\mathrm{0}}=z_{\mathrm{0},\mathrm{v}}/z_{\mathrm{0},\mathrm{vh}}$, which is a SURFEX input parameter that is used to diagnose the roughness for heat and vapor fluxes from the prescribed momentum roughness length, z_0,v, is tested. The lower r_z0 is, the higher the turbulent fluxes become; in ISBA, the default value is r_z0=10, while it is r_z0= exp(1) in MEB (Napoly et al., 2017) following Lo (1995) and Yang and Friedl (2003), who propose values more adapted for forest covers. The uncertainty associated with this parameter motivated this sensitivity test, and values from 1 to 10 were tested. Note that for certain local-scale studies with ISBA, values in excess of 10 have been used; however, for current applications in hydrology and atmospheric modeling in SURFEX, the default of 10 is used so this is the limit used herein. As it turns out, the sensitivity to this parameter during the snow period is very low for the three studied sites, with a maximum variation of 3 W in the RMSE in H and LW↑ compared to the results obtained using the default values. Therefore, no modification was made.

4.2.4 Snow interception

The snow interception parameterization in MEB is based on Hedstrom and Pomeroy (1998) and the implementation in MEB is described in detail in Boone et al. (2017). This scheme is widely used in LSMs; however, note that some alternatives exist which are significantly different (e.g. Roesch et al., 2001). Therefore in the current study, we investigate the sensitivity of the results to this physical process using a rather radical test for which the simulations are repeated using MEB with the maximum interception storage set to zero, thereby effectively turning “off” the snow interception and loading parameterization. It was found that for these particular sites, the process of snow interception by the canopy vegetation has only a mild impact on the snowpack below the canopy and a fairly small impact on the fluxes to the atmosphere. The impact of removing the snow interception on simulated 12-year average annual cycle of snow depth is shown in Fig. 5 in which the MEB simulation without this process is represented by the dashed red curve. The RMSE of the simulated snow depth varies by less than 10 % between this test and the default MEB simulation averaged over the three sites. The maximum snow peak is increased on average by 10 % (3 cm), and the total LE is decreased by approximately 4 % on average for all three sites resulting in small compensating increases in H. Perhaps most significantly, there is virtually no impact on the last day of snow score defined in Sect. 4.1.2, which corresponds to the main difference (and improvement) when comparing MEB to ISBA. We conclude that even if snow interception is a key physical process in land surface schemes such as MEB (Rutter et al., 2009), it does not have much of an impact for the BERMS sites in terms of the improvement of the snowpack modeling, although it does have a relatively large impact on the sublimation, but no specific observations of this flux are available for the sites herein. Therefore, we leave the parameterization with its default parameter values pending the results of future studies.

Table 8 The RMSE and BIAS for the ISBA-MEB and ISBA experiments for snow depth calculated using half-hourly values for which snow is present in the measurements.

5 Conclusions

The impact of snow conditions on the surface fluxes and state variables simulated using the multi-energy balance (MEB) option, which has been recently implemented in the ISBA LSM on the SURFEX platform, is evaluated in this study. The default representation of the surface energy balance in ISBA consists of a single composite soil–vegetation layer for which physical parameters for the two surfaces are weighted by a fraction of surface covered by vegetation. The new option improves the representation of forests through the addition of two explicit layers: a bulk vegetation canopy and a forest surface litter layer. A new energy budget is computed for the bulk canopy, while the below-canopy energy budget is computed for the litter layer. The evaluation has been carried out using 12 years of observations available from the three BERMS (Boreal Ecosystem Research Study) experimental sites that have been used in numerous studies (e.g. Bartlett et al., 2006) and the recent ESM-SnowMIP intercomparison study (Krinner et al., 2018; Menard et al., 2020) and can be considered as a benchmark for evaluating LSMs simulating cold season processes for forested areas.

During periods without snow cover, comparable results and conclusions from a previous study (Napoly et al., 2017) that compared ISBA with ISBA-MEB were confirmed. They can be summarized as follows: due to the shading effect of the canopy layer and the low thermal diffusivity of the litter layer, the ground heat flux daily amplitude, as well as soil temperature daily amplitudes, is significantly reduced. The result is that the ground heat flux was found to be in much better agreement with the measurements (RMSE of 47.1 W m⁻² with ISBA versus 10.9 W m⁻² with MEB) over the entire 12-year integration period. The reduced energy used for conductive heat flux (note that net radiation is barely impacted between the two versions) tends to be manifested as concomitant increases in the daily peak sensible heat flux (RMSE of 58.9 W m⁻² with ISBA and 48.4 W m⁻² with MEB) as the latent heat flux is limited by the evaporative demand. In spring, the latent heat flux is also improved mainly owing to more a limited contribution of the ground evaporation due to the addition of a litter layer (main effect), a decreased surface roughness, and lower wind speeds compared to ISBA. Available measurements of shortwave radiation below the canopy also show the ability of MEB to model the radiative transfer through the canopy.

During snow periods, MEB provides an improved realism of the decoupling between the atmosphere and the ground below the snowpack. Since MEB has eliminated the fractional burying of the ground-based snowpack by the vegetation layer (the p_snv parameterization), the below-canopy surface relies on p_sng uniquely and therefore the ground can be completely covered and insulated for a relatively shallow total snow accumulation (the default value is 0.10 m, which was used in the current study). With the default composite version of ISBA, this decoupling cannot be represented as an effective snow fraction and is calculated in the range of 10 % to 20 %. Consequently, a large fraction of the surface is directly connected with the atmosphere and ground heat flux becomes large (negative directed downward) when the atmospheric temperature decreases. This leads to a strong unrealistic cooling of the soil with an average bias for the all sites, depths and years of −5 compared to −0.1 ^∘C with MEB. In addition to the insulating effect of the litter, this improved representation of the ground heat flux provides energy to turbulent fluxes (and again mostly the sensible heat flux) which were underestimated. The average RMSE for the sensible heat flux calculated for snow periods drops from 54.9 to 45.0 W m⁻².

The improved soil temperature simulation (reduced cold bias) was also found throughout the soil column, and was confirmed by the observations that extended to a 1 m soil depth. It should be noted that this significantly impacts soil phase changes in the model, and thus the modeling of the permafrost in forest regions within large-scale coupled-atmospheric or hydrological simulations, for example. This aspect is the subject of current research and is beyond the scope of this study, but here we note that ISBA coupled to a river routing scheme tends to simulate a river discharge peak, owing to springtime thaw and snowmelt in historical offline simulations north of 50^∘ N, that is too early in comparison to observations and this is attributed to a precocious snow melt and soil thaw in boreal forest regions Decharme et al. (2019).

The impact of the explicit vegetation canopy and litter layer on the snowpack simulation is significant. In general, the snow depth is improved with MEB. The average RMSE calculated for all sites and years is 5.1 cm for MEB, while with ISBA it is 9.1 cm. This is due to a general better agreement during the entire season and specifically from better representation during the melting period. Indeed, the snow melts, on average, 24 d too early with ISBA, while the melt occurs only 1 d early with MEB. This mainly arises in ISBA due to the more direct coupling of the ground with the atmosphere. When the air temperature increases above freezing, the composite-layer temperature also warms, thereby heating the snowpack from below and provoking melt. In MEB, the snowpack occupies the whole fraction of ground so that it can only melt from its surface due to a positive energy balance excess. In addition, because of the lower fractional coverage of snow in ISBA, sublimation represents only 2 % of the total snowfall loss, while it represents approximately 27 % when using MEB. About half of this quantity corresponds to snow intercepted by the canopy, the other half is directly sublimated from the snowpack. While there are no direct estimations of sublimation available at the BERMS sites, it is found that these values correspond well with the total sublimation and partitioning values quoted in the literature for forested sites.

Two hydrological impacts are to be expected with the MEB option. First, the soil does not freeze as deep or for as long. With ISBA, the deep soil can be frozen to well below 1 m depth more than half of the year, while the depth of the 0 ^∘C isotherm and the soil water frozen fraction are both considerably less with MEB. At a soil depth of 1 m, on average over all sites and years, the daily temperature falls below 0 ^∘C on 33 d in the observations, while it is simulated as 35 d per year with MEB compared to 188 d with ISBA. This effect tends to cause ISBA to have a later peak in total runoff. It will be studied in global runs coupled with a hydrological model in the near future.

It should be noted that when doing local-scale uncoupled studies with ISBA, the default p_snv parameterization can be changed such that it rapidly reaches unity after a few centimeters of snow cover has developed. The result is a simulation that is more consistent with MEB in terms of soil temperature and snowpack duration. However, such a configuration can not be used in coupled (with an atmospheric model) mode since there would be a huge positive bias in the upwelling shortwave radiation, notably in spring, thereby potentially having an impact on the simulated high-latitude radiative feedback (which is important to simulate correctly for climate change prediction). In contrast, if the standard p_snv parameterization is used in ISBA (similar to the coupled configuration), ISBA tends to produce the effects cited herein. There, MEB has removed this inconsistency.

For application of MEB in spatially distributed applications, such as for numerical weather prediction (NWP) or hydrological forecasting, such parameters are referred to as “primary parameters” and are generally fixed or prescribed from look-up tables based upon land use classification or plant functional types (set in the SURFEX physiographic database ECOCLIMAP): thus, care must be taken to define values and explore model sensitivity. This is in contrast to so-called “secondary parameters”, which are derived based upon input primary parameters or input physiographic data (such as LAI, soil texture, etc.). As it turns out, the model showed significant sensitivity to only one parameter among the three tested, which was the longwave radiation transmission coefficient. In this study, a slightly lower value of 0.4 was finally chosen compared to the default value of 0.5 from (Boone et al., 2017). The default values of the other parameters from the aforementioned study were unchanged as a result of this study.

Work is currently underway to use MEB within ISBA for forest covers in many of the applications using SURFEX. For regional (covering France) high-resolution (8 km grid) hydrological and surface state forecasting and analysis, Météo-France uses the SURFEX-ISBA-MODCOU hydrometeorological model version 2 (SIM2: Le Moigne et al., 2020), and MEB is being tested for future implementation. There are also preparations underway to use MEB in the coupled SURFEX–CTRIP system in both offline mode and coupled to CNRM-CM (Decharme et al., 2019). Work is also underway at the recently instrumented Col de Porte forest site (Lejeune et al., 2019; Helbig et al., 2020) to use MEB coupled to the detailed snow process model CROCUS (Vionnet et al., 2012), which is used for, among many applications and fundamental research, operational avalanche prediction for mountainous areas in France. This is also a way to test the model in the relatively warm and wet climate that is typical of the French Alps and very different from the BERMS sites. There are longer-term plans to use MEB in the operational regional and global meteorological prediction models AROME and ARPEGE, respectively. Thus, we will continue to evaluate MEB from the local scale (in order to study processes in detail) up to global scales in both offline and coupled land–atmosphere–hydrology model platforms.