A Flexible Snow Model (FSM 2.1.1) including a forest canopy

Multiple options for representing physical processes in forest canopies are added to FSM, which is a model with multiple options for representing physical processes in snow on the ground. The canopy processes represented are shortwave and longwave radiative transfer; turbulent transfers of heat and moisture; and interception, sublimation, unloading, and melt of snow in the canopy. There are options for Beer's law or two-stream approximation canopy radiative transfer, linear or non-linear canopy snow interception efficiency, and time- and melt-dependent or temperature- and wind-dependent canopy snow unloading. Canopy mass and energy balance equations can be solved with one or two model layers. Model behaviour on stand scales is compared with observations of above- and below-canopy shortwave and longwave radiation, below-canopy wind speed, snow mass on the ground, and subjective estimates of canopy snow load. Large-scale simulations of snow cover extent, snow mass, and albedo for the Northern Hemisphere are compared with observations and land-only simulations by state-of-the-art Earth system models. Without accounting for uncertainty in forest structure metrics and parameter values, the ranges of multi-physics ensemble simulations are not as wide as seen in intercomparisons of existing models. FSM2 provides a platform for rapid investigation of sensitivity to model structure and parameter values or ensemble-based data assimilation for snow in open and forested environments.

1 Introduction

Several snow models based on physical principles of energy and mass conservation but offering a range of alternative parametrizations for uncertain energy and mass exchange processes have been developed recently for cryospheric and hydrological applications (Clark et al., 2015; Lafaysse et al., 2017; Niu et al., 2011; Sauter et al., 2020). The Factorial Snow Model (FSM; Essery, 2015) is one such multi-physics model that has proved popular with users because the code is freely available, compact, easy to use, flexible, and thoroughly documented. All of the model state variables and parameters are made accessible through restart and control files for model calibration and data assimilation. At the time of writing, FSM has been used in more than 30 peer-reviewed publications and seven postgraduate theses by users in 11 countries (Essery, 2024). Applications have included snow data assimilation (Alonso-González et al., 2022); evaluation of snow simulations in the European Alps (Magnusson et al., 2015), Norway (Magnusson et al., 2019), and the western Himalaya (Pritchard, 2020); and construction of gridded snow datasets for the Iberian Peninsula (Alonso-González et al., 2018).

FSM was originally intended for investigating the range of results from existing physically based models in simulating snow accumulation and melt on land. Parametrizations of five important processes (decreasing albedo and increasing density of snow with age, increasing thermal conductivity with snow density, storage of liquid water in snow, and suppression of turbulent fluxes in stably stratified atmospheric surface layers) can be switched on or off independently, giving up to 32-member ensembles of simulations. Although neglecting any of these processes is theoretically expected to give poor results, options to neglect them were included because they are neglected in some existing snow models. In fact, Günther et al. (2020) have subsequently shown that it can be impossible to distinguish between model configurations that include or neglect specific processes in the face of parameter uncertainty and with limited evaluation data.

The original version of FSM allowed for partial snow cover on the ground but did not account for exposed vegetation above snow. Large areas of the Northern Hemisphere have both forests and seasonal snow cover, and influences of forest–snow interactions on weather and climate have been of long-standing interest (Chalita and Le Treut, 1994; Thomas and Rowntree, 1992; Viterbo and Betts, 1999). Forest canopies can intercept falling snow, shade underlying snow surfaces from solar radiation but increase incoming thermal radiation, increase drag on the atmosphere and reduce turbulent exchanges between underlying snow and the air, drop litter that decreases the albedo of underlying snow, and mask the albedo of snow-covered land. Intercepted snow may unload to the ground, melt in the canopy, or sublimate back to the atmosphere, but Lundquist et al. (2021) found large differences between results from common parametrizations of these processes that are based on limited observations. Qu and Hall (2014) found large differences in snow–albedo feedbacks for global climate models with differing representations of snow albedo masking by forests. Climate models have to use simple parametrizations for computational efficiency; they often represent forests as a single homogeneous layer or simply modify parameters describing the aerodynamic and radiative properties of the land surface according to vegetation cover. Because a single model layer cannot represent vertical temperature gradients and decoupling between solar radiation absorption in the upper canopy and thermal emission from the lower canopy, there has been recent interest in canopy models using two layers (Gouttevin et al., 2015; Todt et al., 2018) or more (Ryder et al., 2016; Bonan et al., 2018). This is a revival of much earlier work on canopy modelling (e.g. Kondo and Watanabe, 1992; Yamazaki et al., 1992) that did not immediately translate to global climate models, which commonly have multi-layer soil and snow models but still use single-layer (big leaf) canopy models (Bonan et al., 2021).

This paper describes and demonstrates options added to FSM for representing interactions between snow and forest canopies. To retain the acronym but to emphasize flexibility over the capability for factorial experiment designs with every possible combination of options, this version 2 has been renamed as the Flexible Snow Model (FSM2). The multi-layer snow model used by both FSM and FSM2 is described in Essery (2015). Although many snow models exist, Essery et al. (2012) noted that a few different process parametrizations are used time and again in different combinations in many of these models; the same can be said of forest canopy models (Lundquist et al., 2021). As in FSM, FSM2 allows switching between parametrization options to improve understanding of how they operate together in a complete energy and mass balance snow model. FSM2 can be run with forest canopies represented by one or two model layers and simple or more sophisticated parametrizations of canopy radiative transfer, snow interception, and unloading. Large-scale climate model land surface schemes with vegetation canopy representations of less or similar complexity to FSM2 include CLASS (Verseghy et al., 1993), CLM (Lawrence et al., 2019), the multi-energy balance (MEB) component of ISBA (Boone et al., 2017), MATSIRO (Takata et al., 2003), the MOSES canopy model (Essery et al., 2003) implemented in JULES (Best et al., 2011), SiB (Sellers et al., 1986), and VISA (Niu and Yang, 2004), but FSM2 does not represent canopy photosynthesis and evapotranspiration. Distributions of snow in forests on small scales (<10 m) have been investigated using SnowPALM (Broxton et al., 2015) and an early release FSM 2.0.1 (Essery, 2019; Mazzotti et al., 2020a). After describing substantial developments in FSM2 since version 2.0.1 in Sect. 2 of this paper, the influences of canopy model options on snow simulations at stand and hemispheric scales are demonstrated with examples in Sect. 3. Results are discussed in relation to other studies in Sect. 4, and this paper concludes with an outlook on opportunities for applications and developments of FSM2.

2 Model description

2.1 Forest and canopy model structure

Bulk forest structure is defined in FSM2 by canopy height h_c, canopy base height h_b, and effective vegetation area index (VAI) Λ including leaves and stems (models that represent transpiration or vegetation dynamics use separate leaf and stem area indices, but FSM2 combines them). Vegetation density is assumed to be constant with height between the base and the top of the canopy. For a two-layer canopy model, the fraction of the canopy in the upper layer is set by parameter f_Λ (0.5 by default), so the midpoints of the upper and lower layers are at heights

$$\begin{array}{}\text{(1)}& z_{\mathrm{1}}=h_{\mathrm{b}}+\left(\mathrm{1}-{\displaystyle \frac{f_{\mathrm{\Lambda }}}{\mathrm{2}}}\right)(h_{\mathrm{c}}-h_{\mathrm{b}})\end{array}$$

and

$$\begin{array}{}\text{(2)}& z_{\mathrm{2}}=h_{\mathrm{b}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{1}-f_{\mathrm{\Lambda }})(h_{\mathrm{c}}-h_{\mathrm{b}}),\end{array}$$

and the layers have vegetation area indices Λ₁=f_ΛΛ and ${\mathrm{\Lambda }}_{\mathrm{2}}=(\mathrm{1}-f_{\mathrm{\Lambda }})\mathrm{\Lambda }$. Rather than being a physically meaningful and species-dependent canopy base height, h_b is an effective height (2 m by default) for a transition from exponential to logarithmic wind speed profiles below the canopy.

The heat capacity of vegetation depends on biomass and water content, but models vary widely in how they determine heat capacity from canopy characteristics, as illustrated in Fig. 1. CLM and MATSIRO neglect canopy heat capacity. SiB and MEB have a low default heat capacity for dry vegetation that is greatly increased by including the heat capacity of intercepted water. MOSES and Gouttevin et al. (2015) have higher heat capacities including separate contributions of leaves and trunks. VISA has an even higher canopy heat capacity, parametrized as a linear function of combined leaf and stem area indices and the mass of intercepted water. In FSM2, the snow-free vegetation canopy heat capacity is C_v=C_ΛΛ, with the default parameter value $C_{\mathrm{\Lambda }}=\mathrm{3.6}\times {\mathrm{10}}^{\mathrm{4}}$ J K⁻¹ m⁻² chosen to give similar values to Gouttevin et al. (2015). The heat capacity of intercepted snow is added for a combined calculation of the energy balance of vegetation and snow in the canopy.

Figure 1 Heat capacity of snow-free canopies as functions of vegetation area index in MEB, MOSES, VISA, and FSM2. Dots show heat capacities calculated using the method of Gouttevin et al. (2015) with leaf and stem area indices for study sites collated by Todt et al. (2018).

2.2 Shortwave radiative transfer

Shortwave radiative transfer in canopies has often been represented as an infinite sum of transmissions and reflections between a single canopy layer and the ground (e.g. Blyth et al., 1999; Stähli et al., 2009; Tribbeck et al., 2004), but a matrix formulation is much easier to generalize to multi-layer models (Zhao and Qualls, 2005). Nomenclature for downward and upward shortwave radiation fluxes at layer boundaries for one-layer and two-layer canopy models is shown in Fig. 2. Each layer has reflectivities R_b and R_d and transmissivities τ_b and τ_d for direct-beam and diffuse radiation, respectively, and forward-scattering fraction s_b for direct-beam radiation. Forward scattering and reflections from the canopy and the snow or ground surface (albedo α) are assumed to be diffuse. Depending on the availability of measurements for driving FSM2, diffuse and direct-beam shortwave radiation components S_↓dif and S_↓dir above the canopy are read as inputs or global radiation is divided into components using the method of Erbs et al. (1982); this parametrizes the diffuse fraction as an empirical function of the ratio between global radiation at the surface and the top of the atmosphere. The optical properties of canopy layers can either be set by bulk parameters in a Beer's law option or calculated from the properties of individual canopy elements in a two-stream approximation. Random orientations are assumed for the canopy elements in either option, although this could be generalized (Otto and Trautmann, 2008). Beer's law and the two-stream approximation can give similar predictions of broadband canopy albedo and transmission for surface energy balance calculations (Essery, 2013), but the two-stream approximation is more accurate for spectrally resolved calculations (Wang, 2003).

Figure 2 Shortwave radiation fluxes, sensible heat fluxes, temperatures, and aerodynamic resistances in one-layer (a) and two-layer (b) canopy models. Arrows show the directions in which fluxes are defined to be positive. Longwave radiation and moisture fluxes are numbered and directed in the same ways as shortwave radiation and sensible heat fluxes, respectively.

2.2.1 One-layer canopy model

Downward and upward diffuse shortwave radiation fluxes at the top and bottom of a single canopy layer are related by

$$\begin{array}{}\text{(3)}& {\displaystyle }{S}_{\downarrow \mathrm{1}}& {\displaystyle }=R_{\mathrm{d}}{S}_{\uparrow \mathrm{1}}+{\mathit{\tau }}_{\mathrm{d}}{S}_{\downarrow \mathrm{dif}}+s_{\mathrm{b}}{S}_{\downarrow \mathrm{dir}},\text{(4)}& {\displaystyle }{S}_{\uparrow \mathrm{1}}& {\displaystyle }=\mathit{\alpha }{S}_{\downarrow \mathrm{1}}+\mathit{\alpha }{\mathit{\tau }}_{\mathrm{b}}{S}_{\downarrow \mathrm{dir}},\text{(5)}& {\displaystyle }{S}_{\uparrow \mathrm{0}}& {\displaystyle }={\mathit{\tau }}_{\mathrm{d}}{S}_{\uparrow \mathrm{1}}+R_{\mathrm{d}}{S}_{\downarrow \mathrm{dif}}+R_{\mathrm{b}}{S}_{\downarrow \mathrm{dir}},\end{array}$$

which can be written as a matrix equation

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& \left(\begin{array}{ccc}\mathrm{1}& -R_{\mathrm{d}}& \mathrm{0}\\ -\mathit{\alpha }& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& -{\mathit{\tau }}_{\mathrm{d}}& \mathrm{1}\end{array}\right)\left(\begin{array}{c}{S}_{\downarrow \mathrm{1}}\\ S_{\uparrow \mathrm{1}}\\ S_{\uparrow \mathrm{0}}\end{array}\right)=\\ & \left(\begin{array}{c}{\mathit{\tau }}_{\mathrm{d}}\\ \mathrm{0}\\ R_{\mathrm{d}}\end{array}\right)S_{\downarrow \mathrm{dif}}+\left(\begin{array}{c}{s}_{\mathrm{b}}\\ \mathit{\alpha }{\mathit{\tau }}_{\mathrm{b}}\\ R_{\mathrm{b}}\end{array}\right)S_{\downarrow \mathrm{dir}}.\end{array}\end{array}$$

With fluxes obtained by solving this equation, net shortwave radiation,

$$\begin{array}{}\text{(7)}& {\mathrm{SW}}_{\mathrm{v}}=S_{\downarrow \mathrm{dif}}-S_{\downarrow \mathrm{1}}+S_{\uparrow \mathrm{1}}-S_{\uparrow \mathrm{0}}+(\mathrm{1}-{\mathit{\tau }}_{\mathrm{b}})S_{\downarrow \mathrm{dir}},\end{array}$$

is absorbed by the vegetation and

$$\begin{array}{}\text{(8)}& {\mathrm{SW}}_{\mathrm{s}}=(\mathrm{1}-\mathit{\alpha })(S_{\downarrow \mathrm{1}}+{\mathit{\tau }}_{\mathrm{b}}{S}_{\downarrow \mathrm{dir}})\end{array}$$

is absorbed by the snow or ground surface.

2.2.2 Two-layer canopy model

The matrix equation for diffuse shortwave radiation fluxes at the boundaries of two canopy layers is

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& \left(\begin{array}{ccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& -R_{\mathrm{d},\mathrm{1}}& \mathrm{0}\\ -{\mathit{\tau }}_{\mathrm{d},\mathrm{2}}& \mathrm{1}& -R_{\mathrm{d},\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\mathit{\alpha }& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ -R_{\mathrm{d},\mathrm{2}}& \mathrm{0}& -{\mathit{\tau }}_{\mathrm{d},\mathrm{2}}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& -{\mathit{\tau }}_{\mathrm{d},\mathrm{1}}& \mathrm{1}\end{array}\right)\left(\begin{array}{c}{S}_{\downarrow \mathrm{1}}\\ S_{\downarrow \mathrm{2}}\\ S_{\uparrow \mathrm{2}}\\ S_{\uparrow \mathrm{1}}\\ S_{\uparrow \mathrm{0}}\end{array}\right)=\\ & \left(\begin{array}{c}{\mathit{\tau }}_{\mathrm{d},\mathrm{1}}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ R_{\mathrm{d},\mathrm{1}}\end{array}\right)S_{\downarrow \mathrm{dif}}+\left(\begin{array}{c}{s}_{\mathrm{b},\mathrm{1}}\\ s_{\mathrm{b},\mathrm{2}}{\mathit{\tau }}_{\mathrm{b},\mathrm{1}}\\ \mathit{\alpha }{\mathit{\tau }}_{\mathrm{b},\mathrm{1}}{\mathit{\tau }}_{\mathrm{b},\mathrm{2}}\\ R_{\mathrm{b},\mathrm{2}}{\mathit{\tau }}_{\mathrm{b},\mathrm{1}}\\ R_{\mathrm{b},\mathrm{1}}\end{array}\right)S_{\downarrow \mathrm{dir}}.\end{array}\end{array}$$

Net shortwave radiation absorbed by vegetation in the two layers and by the snow or ground surface is

$$\begin{array}{}\text{(10)}& {\displaystyle }{\mathrm{SW}}_{\mathrm{v},\mathrm{1}}=S_{\downarrow \mathrm{dif}}-S_{\downarrow \mathrm{1}}+S_{\uparrow \mathrm{1}}-S_{\uparrow \mathrm{0}}+(\mathrm{1}-{\mathit{\tau }}_{\mathrm{b},\mathrm{1}})S_{\downarrow \mathrm{dir}},\text{(11)}& {\displaystyle }{\mathrm{SW}}_{\mathrm{v},\mathrm{2}}=S_{\downarrow \mathrm{1}}-S_{\downarrow \mathrm{2}}+S_{\uparrow \mathrm{2}}-S_{\uparrow \mathrm{1}}+{\mathit{\tau }}_{\mathrm{b},\mathrm{1}}(\mathrm{1}-{\mathit{\tau }}_{\mathrm{b},\mathrm{2}})S_{\downarrow \mathrm{dir}},\end{array}$$

and

$$\begin{array}{}\text{(12)}& {\mathrm{SW}}_{\mathrm{s}}=(\mathrm{1}-\mathit{\alpha })(S_{\downarrow \mathrm{2}}+{\mathit{\tau }}_{\mathrm{b},\mathrm{1}}{\mathit{\tau }}_{\mathrm{b},\mathrm{2}}{S}_{\downarrow \mathrm{dir}}).\end{array}$$

2.2.3 Beer's law option

The fraction of radiation incident from above at elevation angle θ transmitted without interception through canopy layer n is parametrized as

$$\begin{array}{}\text{(13)}& {\mathit{\tau }}_{\mathrm{b},n}=\exp (-k_{\mathrm{ext}}{\mathrm{\Lambda }}_n/\sin \mathit{\theta }),\end{array}$$

with extinction coefficient k_ext=0.5 by default for randomly oriented canopy elements. Integrating Eq. (13) over the sky hemisphere to find the transmission of diffuse radiation through a layer without interception results in an exponential integral (Nijssen and Lettenmaier, 1999) that can be closely approximated by

$$\begin{array}{}\text{(14)}& {\mathit{\tau }}_{\mathrm{d},n}=\exp (-\mathrm{1.6}{k}_{\mathrm{ext}}{\mathrm{\Lambda }}_n).\end{array}$$

Transmission is thus higher for diffuse than direct-beam radiation for elevation angles less than 39°. Forward scattering is neglected ($s_{\mathrm{b},n}=\mathrm{0}$). Diffuse and direct-beam canopy layer reflectivities are taken to be $R_{\mathrm{d},n}=(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d},n}){\mathit{\alpha }}_{\mathrm{c}}$ and $R_{\mathrm{b},n}=(\mathrm{1}-{\mathit{\tau }}_{\mathrm{b},n}){\mathit{\alpha }}_{\mathrm{c}}$ for dense-canopy albedo α_c in the limit Λ→∞. For a canopy layer with snow cover fraction f_cs, this is given by

$$\begin{array}{}\text{(15)}& {\mathit{\alpha }}_{\mathrm{c}}=(\mathrm{1}-f_{\mathrm{cs}}){\mathit{\alpha }}_{c\mathrm{0}}+f_{\mathrm{cs}}{\mathit{\alpha }}_{\mathrm{cs}}\end{array}$$

for snow-free and snow-covered dense-canopy albedo parameters α_c0 and α_cs, with default values 0.1 and 0.3, respectively, based on Bartlett and Verseghy (2015).

2.2.4 Two-stream approximation option

From a review of two-stream radiative transfer approximations by Meador and Weaver (1980), Dickinson (1983) adapted the hemispheric constant method for isotropic multiple scattering of light in a homogeneous canopy layer. This was used by Sellers et al. (1986) in SiB and is now used in CLM.

The two-stream model can account for transmission of light through leaves, but FSM2 assumes that canopy elements are opaque and have reflectivities α_Λ0 when snow-free and α_Λs when fully covered with snow. Partial canopy snow cover gives

$$\begin{array}{}\text{(16)}& {\mathit{\alpha }}_{\mathrm{\Lambda }}=(\mathrm{1}-f_{\mathrm{cs}}){\mathit{\alpha }}_{\mathrm{\Lambda }\mathrm{0}}+f_{\mathrm{cs}}{\mathit{\alpha }}_{\mathrm{\Lambda }\mathrm{s}}.\end{array}$$

Solutions of the two-stream equations from Meador and Weaver (1980) involve coefficients

$$\begin{array}{}\text{(17)}& {\mathit{\gamma }}_{\mathrm{1}}=\mathrm{2}[\mathrm{1}-(\mathrm{1}-\mathit{\beta }\left)\mathit{\omega }\right],{\mathit{\gamma }}_{\mathrm{2}}=\mathrm{2}\mathit{\beta }\mathit{\omega },{\mathit{\gamma }}_{\mathrm{3}}={\mathit{\beta }}_{\mathrm{0}},{\mathit{\gamma }}_{\mathrm{4}}=\mathrm{1}-{\mathit{\beta }}_{\mathrm{0}}.\end{array}$$

For flat, opaque, and randomly oriented leaves, ω=α_Λ is the fraction of incident radiation that is scattered, $\mathit{\beta }=\mathrm{2}/\mathrm{3}$ is the fraction of scattered diffuse radiation that is directed back into the upward hemisphere, and

$$\begin{array}{}\text{(18)}& {\mathit{\beta }}_{\mathrm{0}}=(\mathrm{0.5}+\mathit{\mu })\left[\mathrm{1}-\mathit{\mu }\ln \left({\displaystyle \frac{\mathrm{1}+\mathit{\mu }}{\mathit{\mu }}}\right)\right]\end{array}$$

is the upscatter fraction for direct-beam radiation with μ=sin θ. The reflectivity and transmissivity of layer n for diffuse radiation are

$$\begin{array}{}\text{(19)}& R_{\mathrm{d},n}={\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{2}}(\mathrm{1}-e^{-\mathrm{2}kl})}{k+{\mathit{\gamma }}_{\mathrm{1}}+(k-{\mathit{\gamma }}_{\mathrm{1}})e^{-\mathrm{2}kl}}}\end{array}$$

and

$$\begin{array}{}\text{(20)}& {\mathit{\tau }}_{\mathrm{d},n}={\displaystyle \frac{\mathrm{2}k{e}^{-kl}}{k+{\mathit{\gamma }}_{\mathrm{1}}+(k-{\mathit{\gamma }}_{\mathrm{1}})e^{-\mathrm{2}kl}}}\end{array}$$

for extinction coefficient $k=({\mathit{\gamma }}_{\mathrm{1}}^{\mathrm{2}}-{\mathit{\gamma }}_{\mathrm{2}}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{2}}$ and optical thickness l=k_extΛ_n. The direct-beam reflectivity, forward-scattering fraction, and transmissivity are

$$\begin{array}{}\text{(21)}& {\displaystyle }{R}_{\mathrm{b},n}={\displaystyle \frac{\begin{array}{c}\mathit{\omega }\left[\right(\mathrm{1}-k\mathit{\mu }\left)\right({\mathit{\alpha }}_{\mathrm{2}}+k{\mathit{\gamma }}_{\mathrm{3}})e^{kl}-(\mathrm{1}+k\mathit{\mu })\\ ({\mathit{\alpha }}_{\mathrm{2}}-k{\mathit{\gamma }}_{\mathrm{3}})e^{-kl}-\mathrm{2}k({\mathit{\gamma }}_{\mathrm{3}}-{\mathit{\alpha }}_{\mathrm{2}}\mathit{\mu })e^{-l/\mathit{\mu }}]\end{array}}{(\mathrm{1}-k^{\mathrm{2}}{\mathit{\mu }}^{\mathrm{2}})\left[\right(k+{\mathit{\gamma }}_{\mathrm{1}})e^{kl}+(k-{\mathit{\gamma }}_{\mathrm{1}}\left)e^{-kl}\right]}},\text{(22)}& {\displaystyle }{s}_{\mathrm{b},n}={\displaystyle \frac{\begin{array}{c}\mathit{\omega }{e}^{-l/\mathit{\mu }}\left[\right(\mathrm{1}-k\mathit{\mu }\left)\right({\mathit{\alpha }}_{\mathrm{1}}-k{\mathit{\gamma }}_{\mathrm{4}})e^{-kl}-(\mathrm{1}+k\mathit{\mu })\\ ({\mathit{\alpha }}_{\mathrm{1}}+k{\mathit{\gamma }}_{\mathrm{4}})e^{kl}]+\mathrm{2}k\mathit{\omega }({\mathit{\gamma }}_{\mathrm{4}}+{\mathit{\alpha }}_{\mathrm{1}}\mathit{\mu })\end{array}}{(\mathrm{1}-k^{\mathrm{2}}{\mathit{\mu }}^{\mathrm{2}})\left[\right(k+{\mathit{\gamma }}_{\mathrm{1}})e^{kl}+(k-{\mathit{\gamma }}_{\mathrm{1}}\left)e^{-kl}\right]}},\end{array}$$

and

$$\begin{array}{}\text{(23)}& {\mathit{\tau }}_{\mathrm{b},n}=e^{-l/\mathit{\mu }},\end{array}$$

where ${\mathit{\alpha }}_{\mathrm{1}}={\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{4}}+{\mathit{\gamma }}_{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{3}}$ and ${\mathit{\alpha }}_{\mathrm{2}}={\mathit{\gamma }}_{\mathrm{1}}{\mathit{\gamma }}_{\mathrm{3}}+{\mathit{\gamma }}_{\mathrm{2}}{\mathit{\gamma }}_{\mathrm{4}}$.

The albedos of snow and vegetation differ strongly between visible and near-infrared wavelengths. Spectrally resolved shortwave radiation fluxes are available to land surface schemes when coupled to atmospheric models, but they are rarely available from measurements. FSM2 currently only calculates canopy reflectivities and transmissivities from average broadband albedos, but averages of separate visible and near-infrared calculations differ by less than 0.03 for all combinations of f_cs, θ, and Λ.

Equation (19) gives the dense-canopy albedo for diffuse radiation as

$$\begin{array}{}\text{(24)}& {\mathit{\alpha }}_{\mathrm{c}}={\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{2}}}{k+{\mathit{\gamma }}_{\mathrm{1}}}},\end{array}$$

which is used to select default parameter values α_Λ0=0.27 and α_Λs=0.65 that match the default snow-free and snow-covered dense-canopy albedos used with Beer's law. Because of absorption of multiply reflected light, the canopy albedo is much lower than the reflectivity of individual canopy elements even when they are covered with snow. Canopy transmission, calculated by dividing above-canopy by sub-canopy radiation, is higher than canopy layer transmissivities because of downward reflections. Figure 3 compares canopy albedos and transmission calculated using Beer's law and the two-stream approximation. The most obvious difference is the systematically greater direct-beam transmission calculated with the two-stream approximation. Transmission of diffuse radiation through a snow-free canopy is slightly lower for the two-stream approximation but can be higher when the canopy is snow-covered. Most of the scatter in Fig. 3 comes from differing fractions of canopy snow cover.

Figure 3 Canopy albedo and transmission of diffuse and direct-beam shortwave radiation calculated using Beer's law and the two-stream approximation in 10 000 simulations with randomly selected values of α (0.1–0.8), θ (5–85°), f_cs (0–1), and Λ (0–10). Green points are from simulations with snow-free canopies.

2.3 Longwave radiation

Vegetation layer n and snow or ground surface temperatures are T_v,n and T_s. Transmission of longwave radiation from the atmosphere through the canopy is given by Eq. (14), in the same way as transmission of diffuse solar radiation. Vegetation, snow, and ground emissivities are assumed to be equal to 1, so no longwave radiation is reflected; this makes the radiative transfer equations tractable without resorting to matrix solutions.

2.3.1 One-layer canopy model

Downward and upward longwave radiation fluxes at the top and bottom of a single canopy layer are related by

$$\begin{array}{}\text{(25)}& {\displaystyle }{L}_{\downarrow \mathrm{1}}={\mathit{\tau }}_{\mathrm{d}}{\mathrm{LW}}_{\downarrow }+(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d}})\mathit{\sigma }{T}_{\mathrm{v}}^{\mathrm{4}},\text{(26)}& {\displaystyle }{L}_{\uparrow \mathrm{1}}=\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}},\end{array}$$

and

$$\begin{array}{}\text{(27)}& L_{\uparrow \mathrm{0}}={\mathit{\tau }}_{\mathrm{d}}{L}_{\uparrow \mathrm{1}}+(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d}})\mathit{\sigma }{T}_{\mathrm{v}}^{\mathrm{4}}\end{array}$$

for incoming longwave radiation LW_↓ above the canopy and Stefan–Boltzmann constant σ, from which net longwave radiation

$$\begin{array}{}\text{(28)}& {\mathrm{LW}}_{\mathrm{v}}=(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d}})({\mathrm{LW}}_{\downarrow }-\mathrm{2}\mathit{\sigma }{T}_{\mathrm{v}}^{\mathrm{4}}+\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}})\end{array}$$

is absorbed by the vegetation and

$$\begin{array}{}\text{(29)}& {\mathrm{LW}}_{\mathrm{s}}={\mathit{\tau }}_{\mathrm{d}}{\mathrm{LW}}_{\downarrow }+(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d}})\mathit{\sigma }{T}_{\mathrm{v}}^{\mathrm{4}}-\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}\end{array}$$

is absorbed by the snow or ground surface. Upwelling longwave radiation above the canopy is

$$\begin{array}{}\text{(30)}& L_{\uparrow \mathrm{0}}=(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d}})\mathit{\sigma }{T}_{\mathrm{v}}^{\mathrm{4}}+{\mathit{\tau }}_{\mathrm{d}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}.\end{array}$$

2.3.2 Two-layer canopy model

The longwave radiation absorbed by vegetation in the canopy layers and longwave radiation absorbed by the surface in a two-layer model are

$$\begin{array}{}\text{(31)}& {\displaystyle }\begin{array}{rl}& {\mathrm{LW}}_{\mathrm{v},\mathrm{1}}=(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d},\mathrm{1}})\\ & \left[{\mathrm{LW}}_{\downarrow }-\mathrm{2}\mathit{\sigma }{T}_{\mathrm{v},\mathrm{1}}^{\mathrm{4}}+(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d},\mathrm{2}})\mathit{\sigma }{T}_{\mathrm{v},\mathrm{2}}^{\mathrm{4}}+{\mathit{\tau }}_{\mathrm{d},\mathrm{2}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}\right],\end{array}\text{(32)}& {\displaystyle }\begin{array}{rl}& {\mathrm{LW}}_{\mathrm{v},\mathrm{2}}=(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d},\mathrm{2}})\\ & \left[{\mathit{\tau }}_{\mathrm{d},\mathrm{1}}{\mathrm{LW}}_{\downarrow }+(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d},\mathrm{1}})\mathit{\sigma }{T}_{\mathrm{v},\mathrm{1}}^{\mathrm{4}}-\mathrm{2}\mathit{\sigma }{T}_{\mathrm{v},\mathrm{2}}^{\mathrm{4}}+\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}\right],\end{array}\end{array}$$

and

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{\mathrm{LW}}_{\mathrm{s}}& ={\mathit{\tau }}_{\mathrm{d},\mathrm{1}}{\mathit{\tau }}_{\mathrm{d},\mathrm{2}}{\mathrm{LW}}_{\downarrow }+(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d},\mathrm{1}}){\mathit{\tau }}_{\mathrm{d},\mathrm{2}}\mathit{\sigma }{T}_{\mathrm{v},\mathrm{1}}^{\mathrm{4}}\\ & +(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d},\mathrm{2}})\mathit{\sigma }{T}_{\mathrm{v},\mathrm{2}}^{\mathrm{4}}-\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}.\end{array}\end{array}$$

Upwelling longwave radiation above the canopy is

$$\begin{array}{}\text{(34)}& \begin{array}{rl}{L}_{\uparrow \mathrm{0}}& =(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d},\mathrm{1}})\mathit{\sigma }{T}_{\mathrm{v},\mathrm{1}}^{\mathrm{4}}+(\mathrm{1}-{\mathit{\tau }}_{\mathrm{d},\mathrm{2}}){\mathit{\tau }}_{\mathrm{d},\mathrm{1}}\mathit{\sigma }{T}_{\mathrm{v},\mathrm{2}}^{\mathrm{4}}\\ & +{\mathit{\tau }}_{\mathrm{d},\mathrm{1}}{\mathit{\tau }}_{\mathrm{d},\mathrm{2}}\mathit{\sigma }{T}_{\mathrm{s}}^{\mathrm{4}}.\end{array}\end{array}$$

2.4 Turbulent fluxes

Vertical momentum, sensible heat, and moisture fluxes are parametrized in FSM2 and many other models by integrals of first-order flux–gradient relationships:

$$\begin{array}{}\text{(35)}& {\displaystyle }\mathit{\rho }{u}_{*}^{\mathrm{2}}=\mathit{\rho }{K}_{\mathrm{m}}{\displaystyle \frac{\partial U}{\partial z}},\text{(36)}& {\displaystyle }H=-\mathit{\rho }{c}_p{K}_{\mathrm{H}}{\displaystyle \frac{\partial T}{\partial z}},\end{array}$$

and

$$\begin{array}{}\text{(37)}& E=-\mathit{\rho }{K}_{\mathrm{H}}{\displaystyle \frac{\partial q}{\partial z}},\end{array}$$

where u_* is the friction velocity, ρ and c_p are the density and heat capacity of air, K_m and K_H are eddy diffusivities for momentum and heat or moisture, U is wind speed, T is air temperature, and q is specific humidity. In open areas and above forest canopies (z>h_c), the eddy diffusivities are given by the Prandtl hypothesis,

$$\begin{array}{}\text{(38)}& K_{\mathrm{m},\mathrm{H}}=k{u}_{*}(z-d){\mathit{\varphi }}_{\mathrm{m},\mathrm{H}}^{-\mathrm{1}}\left({\displaystyle \frac{z-d}{L}}\right),\end{array}$$

for von Kármán constant k, displacement height d, Obukhov length L, and similarity functions ϕ_m and ϕ_H described later.

2.4.1 Open areas

Momentum roughness lengths z_0f for snow-free ground and z_0s for snow on fraction f_s of the ground are combined to give a composite surface roughness length:

$$\begin{array}{}\text{(39)}& z_{\mathrm{0}}=z_{\mathrm{0}\mathrm{f}}^{\mathrm{1}-f_{\mathrm{s}}}{z}_{\mathrm{0}\mathrm{s}}^{f_{\mathrm{s}}}.\end{array}$$

The roughness length for heat transfer from the surface is z_0h=0.1z₀, and d=0 in Eq. (38). Integrating Eq. (35) between z₀ (where U=0) and wind measurement height z_U (where U=U_a) in a surface layer with constant u_* gives

$$\begin{array}{}\text{(40)}& u_{*}=k{U}_{\mathrm{a}}{\left[\ln \left({\displaystyle \frac{z_U}{z_{\mathrm{0}}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_U}{L}}\right)+{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{0}}}{L}}\right)\right]}^{-\mathrm{1}},\end{array}$$

where

$$\begin{array}{}\text{(41)}& {\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z}{L}}\right)=\underset{\mathrm{0}}{\overset{z/L}{\int }}\left[{\displaystyle \frac{\mathrm{1}-{\mathit{\varphi }}_{\mathrm{m}}\left(\mathit{\zeta }\right)}{\mathit{\zeta }}}\right]\mathrm{d}\mathit{\zeta }.\end{array}$$

Wind speeds at heights z<z_U are given by

$$\begin{array}{}\text{(42)}& U\left(z\right)={\displaystyle \frac{u_{*}}{k}}\left[\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z}{L}}\right)+{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{0}}}{L}}\right)\right].\end{array}$$

Integrating Eq. (36) between z_0h (where T=T_s) and temperature measurement height z_T (where T=T_a) with constant H gives

$$\begin{array}{}\text{(43)}& H={\displaystyle \frac{\mathit{\rho }{c}_p}{r_{\mathrm{a}}}}(T_{\mathrm{s}}-T_{\mathrm{a}})\end{array}$$

for aerodynamic resistance

$$\begin{array}{}\text{(44)}& r_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{k{u}_{*}}}\left[\ln \left({\displaystyle \frac{z_T}{z_{\mathrm{0}\mathrm{h}}}}\right)-{\mathit{\psi }}_{\mathrm{H}}\left({\displaystyle \frac{z_T}{L}}\right)+{\mathit{\psi }}_{\mathrm{H}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{h}}}{L}}\right)\right]\end{array}$$

and

$$\begin{array}{}\text{(45)}& {\mathit{\psi }}_{\mathrm{H}}\left({\displaystyle \frac{z}{L}}\right)=\underset{\mathrm{0}}{\overset{z/L}{\int }}\left[{\displaystyle \frac{\mathrm{1}-{\mathit{\varphi }}_{\mathrm{H}}\left(\mathit{\zeta }\right)}{\mathit{\zeta }}}\right]\mathrm{d}\mathit{\zeta }.\end{array}$$

Moisture flux,

$$\begin{array}{}\text{(46)}& E={\mathit{\chi }}_{\mathrm{s}}{\displaystyle \frac{\mathit{\rho }}{r_{\mathrm{a}}}}\left[q_{\mathrm{sat}}\right(T_{\mathrm{s}})-q_{\mathrm{a}}],\end{array}$$

is calculated using the same aerodynamic resistance as the sensible heat flux and a moisture availability factor χ_s=1 if E<0 (condensation) or if there is snow on the ground, or

$$\begin{array}{}\text{(47)}& {\mathit{\chi }}_{\mathrm{s}}={\displaystyle \frac{r_{\mathrm{a}}}{r_{\mathrm{a}}+r_{\mathrm{sg}}}}\end{array}$$

otherwise, where r_sg is the resistance for evaporation from soil moisture (a fixed parameter in the absence of interactive soil moisture and photosynthesis models in FSM2). q_a is the specific humidity measured at height z_T, and q_sat(T) is the saturation humidity at temperature T. Latent heat flux LE is calculated by multiplying the moisture flux by the latent heat of sublimation if T_s<T_m or the latent heat of evaporation otherwise.

2.4.2 Forest areas

Aerodynamic resistances and turbulent fluxes above, within, and beneath forest canopies are shown for one-layer and two-layer canopy models in Fig. 2. Resistances r_a, r_v,n, and r_s couple the upper-canopy air space (temperature T_c,1) to the atmosphere, vegetation layers (temperatures T_v,n) to corresponding canopy air space layers, and the ground or snow surface to the lower-canopy air space, respectively. An additional resistance r_c couples the upper and lower-canopy air space layers in the two-layer model. The sensible heat fluxes are parametrized as

$$\begin{array}{}\text{(48)}& H={\displaystyle \frac{\mathit{\rho }{c}_p}{r_{\mathrm{a}}}}(T_{\mathrm{c},\mathrm{1}}-T_{\mathrm{a}})\end{array}$$

from the upper-canopy air space to the atmosphere,

$$\begin{array}{}\text{(49)}& H_{\mathrm{v},n}={\displaystyle \frac{\mathit{\rho }{c}_p}{r_{\mathrm{v},n}}}(T_{\mathrm{v},n}-T_{\mathrm{c},n})\end{array}$$

from vegetation layer n to canopy air space layer n, and

$$\begin{array}{}\text{(50)}& H_{\mathrm{s}}={\displaystyle \frac{\mathit{\rho }{c}_p}{r_{\mathrm{s}}}}(T_{\mathrm{s}}-T_{\mathrm{c},N})\end{array}$$

from the surface to the lower-canopy air space (N=1 for a one-layer canopy model or 2 for a two-layer canopy model). The flux between canopy air space layers in a two-layer model is

$$\begin{array}{}\text{(51)}& H_{\mathrm{c}}={\displaystyle \frac{\mathit{\rho }{c}_p}{r_{\mathrm{c}}}}(T_{\mathrm{c},\mathrm{2}}-T_{\mathrm{c},\mathrm{1}}).\end{array}$$

The corresponding moisture fluxes are

$$\begin{array}{}\text{(52)}& {\displaystyle }E={\displaystyle \frac{\mathit{\rho }}{r_{\mathrm{a}}}}(q_{\mathrm{c},\mathrm{1}}-q_{\mathrm{a}}),\text{(53)}& {\displaystyle }{E}_{\mathrm{v},n}={\mathit{\chi }}_{\mathrm{v},n}{\displaystyle \frac{\mathit{\rho }}{r_{\mathrm{v},n}}}\left[q_{\mathrm{sat}}\right(T_{\mathrm{v},n})-q_{\mathrm{c},n}],\text{(54)}& {\displaystyle }{E}_{\mathrm{s}}={\mathit{\chi }}_{\mathrm{s}}{\displaystyle \frac{\mathit{\rho }}{r_{\mathrm{s}}}}\left[q_{\mathrm{sat}}\right(T_{\mathrm{s}})-q_{\mathrm{c},N}],\end{array}$$

and

$$\begin{array}{}\text{(55)}& E_{\mathrm{c}}={\displaystyle \frac{\mathit{\rho }}{r_{\mathrm{c}}}}(q_{\mathrm{c},\mathrm{2}}-q_{\mathrm{c},\mathrm{1}}),\end{array}$$

where q_c,n is the specific humidity of canopy air space layer n. The moisture availability factor for evaporation from vegetation layer n is ${\mathit{\chi }}_{\mathrm{v},n}=\mathrm{1}$ if $E_{\mathrm{v},n}<\mathrm{0}$ or

$$\begin{array}{}\text{(56)}& {\mathit{\chi }}_{\mathrm{v},n}=f_{\mathrm{cs},n}+(\mathrm{1}-f_{\mathrm{cs},n}){\displaystyle \frac{r_{\mathrm{v},n}}{r_{\mathrm{v},n}+r_{\mathrm{sv}}}}\end{array}$$

otherwise, where r_sv is a fixed resistance for evaporation from snow-free vegetation (100 s m⁻¹ by default, typical of unstressed vegetation).

Displacement height d=0.67h_c and vegetation roughness length z_0v=0.1h_c are used for dense canopies. FSM2 and many other models (e.g. Choudhury and Monteith, 1998; Dolman, 1993; Koivusalo and Kokkonen, 2002; Boone et al., 2017) follow Inoue (1963) in assuming exponential wind profiles within dense canopies. Eddy diffusivities are continuous at the canopy top and have an exponential form,

$$\begin{array}{}\text{(57)}& K\left(z\right)=K\left(h_{\mathrm{c}}\right)\exp \left[\mathit{\eta }\left({\displaystyle \frac{z}{h_{\mathrm{c}}}}-\mathrm{1}\right)\right],\end{array}$$

with η=2.5 by default from Dolman (1993). Within-canopy stability effects are neglected. The wind speed profile above, within, and below the canopy is

$$\begin{array}{}\text{(58)}& \begin{array}{rl}& U\left(z\right)=\\ & \left\{\begin{array}{ll}\frac{u_{*}}{k}\left[\ln \left(\frac{z-d}{z_{\mathrm{0}\mathrm{v}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z-d}{L}\right)+{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{0}\mathrm{v}}}{L}\right)\right]& z\ge h_{\mathrm{c}}\\ U_{\mathrm{c}}\exp \left[\mathit{\eta }\left(\frac{z}{h_{\mathrm{c}}}-\mathrm{1}\right)\right]& h_{\mathrm{b}}<z<h_{\mathrm{c}}\\ U_{\mathrm{b}}\ln \left(\frac{z}{z_{\mathrm{0}}}\right){\left[\ln \left(\frac{h_{\mathrm{b}}}{z_{\mathrm{0}}}\right)\right]}^{-\mathrm{1}}& z_{\mathrm{0}}\le z\le h_{\mathrm{b}},\end{array}\right.\end{array}\end{array}$$

with

$$\begin{array}{}\text{(59)}& \begin{array}{rl}{u}_{*}& =k{U}_{\mathrm{a}}\left[\ln \left({\displaystyle \frac{z_U-d}{z_{\mathrm{0}\mathrm{v}}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_U-d}{L}}\right)\right.\\ & {\left.+{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{v}}}{L}}\right)\right]}^{-\mathrm{1}}.\end{array}\end{array}$$

Continuity is used to calculate wind speeds U_c and U_b at canopy top and base heights h_c and h_b. Integrating Eq. (36) between the relevant heights gives aerodynamic resistances

$$\begin{array}{}\text{(60)}& \begin{array}{rl}{r}_{\mathrm{a}}& ={\displaystyle \frac{\mathrm{1}}{k{u}_{*}}}\left[\ln \left({\displaystyle \frac{z_T-d}{h_{\mathrm{c}}-d}}\right)-{\mathit{\psi }}_{\mathrm{H}}\left({\displaystyle \frac{z_T-d}{L}}\right)\right.\\ & \left.+{\mathit{\psi }}_{\mathrm{H}}\left({\displaystyle \frac{h_{\mathrm{c}}-d}{L}}\right)\right]+{\displaystyle \frac{h_{\mathrm{c}}\left[e^{\mathit{\eta }(\mathrm{1}-z_{\mathrm{1}}/h_{\mathrm{c}})}-\mathrm{1}\right]}{\mathit{\eta }{K}_{\mathrm{H}}\left(h_{\mathrm{c}}\right)}}\end{array}\end{array}$$

for heat transfer between the highest canopy layer and the atmosphere,

$$\begin{array}{}\text{(61)}& r_{\mathrm{c}}={\displaystyle \frac{e^{\mathit{\eta }}{h}_{\mathrm{c}}}{\mathit{\eta }{K}_{\mathrm{H}}\left(h_{\mathrm{c}}\right)}}\left(e^{-\mathit{\eta }{z}_{\mathrm{2}}/h_{\mathrm{c}}}-e^{-\mathit{\eta }{z}_{\mathrm{1}}/h_{\mathrm{c}}}\right)\end{array}$$

between heights z₁ and z₂ within the canopy, and

$$\begin{array}{}\text{(62)}& \begin{array}{rl}{r}_{\mathrm{s}}=& {\displaystyle \frac{\mathrm{1}}{k^{\mathrm{2}}{U}_{\mathrm{b}}}}\ln \left({\displaystyle \frac{h_{\mathrm{b}}}{z_{\mathrm{0}}}}\right)\ln \left({\displaystyle \frac{h_{\mathrm{b}}}{z_{\mathrm{0}\mathrm{h}}}}\right)\\ & +{\displaystyle \frac{e^{\mathit{\eta }}{h}_{\mathrm{c}}}{\mathit{\eta }{K}_{\mathrm{H}}\left(h_{\mathrm{c}}\right)}}\left(e^{-\mathit{\eta }{h}_{\mathrm{b}}/h_{\mathrm{c}}}-e^{-\mathit{\eta }{z}_N/h_{\mathrm{c}}}\right)\end{array}\end{array}$$

between the surface and the lowest canopy layer. The resistance for sensible heat flux from vegetation to the air within canopy layer n is given by

$$\begin{array}{}\text{(63)}& {\displaystyle \frac{\mathrm{1}}{r_{\mathrm{v},n}}}=C_{\mathrm{leaf}}{\mathrm{\Lambda }}_{n}U(z_n{)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

with C_leaf=0.05 m${}^{\mathrm{1}/\mathrm{2}}$ s${}^{-\mathrm{1}/\mathrm{2}}$ by default from Lawrence et al. (2019). Calculating conductance $r_{\mathrm{v}}^{-\mathrm{1}}$ rather than resistance r_v avoids dividing by small numbers as Λ→0. Canopy conductances of this form are stated in several model description papers without reference or simply with reference to earlier models. The $U^{\mathrm{1}/\mathrm{2}}$ dependence is characteristic of engineering expressions for laminar flow over plates and appears in biophysical literature for flow over leaves at least as far back as Raschke (1960).

Wang (2012) pointed out that exponential wind speed profiles do not satisfy no-slip boundary conditions at the ground and do not converge to logarithmic profiles for zero canopy density. Moreover, Inoue (1963) predicted that η should be a function of canopy density, but models generally take it to be a constant parameter. The logarithmic wind speed profile below height h_b in Eq. (58) is commonly adopted in models to impose a no-slip boundary condition. In FSM2, dense-canopy and open conductances are weighted by the vertically projected vegetation fraction $f_{\mathrm{v}}=\mathrm{1}-\exp (-k_{\mathrm{ext}}\mathrm{\Lambda })$ and (1−f_v), respectively, and combined in parallel to get resistances for sparse canopies. Sub-canopy wind speeds are calculated as weighted averages of wind speeds from Eqs. (42) and (58), as shown in Fig. 4.

Figure 4 Wind speeds above, within, and below canopies with vegetation area indices (from right to left) 0, 0.5, 1, 2, 4, and 8.

2.4.3 Stability functions

Atmospheric stability is characterized by Obukhov length

$$\begin{array}{}\text{(64)}& L=-{\displaystyle \frac{\mathit{\rho }{c}_p{T}_{\mathrm{a}}{u}_{*}^{\mathrm{3}}}{kgH}}.\end{array}$$

As in Bonan et al. (2018) (but without modification for a roughness sublayer above canopies), the stability functions used in FSM2 are

$$\begin{array}{}\text{(65)}& {\mathit{\varphi }}_{\mathrm{m}}\left(\mathit{\zeta }\right)=\left\{\begin{array}{ll}(\mathrm{1}-\mathrm{16}\mathit{\zeta }{)}^{-\mathrm{1}/\mathrm{4}}& \mathit{\zeta }<\mathrm{0}\\ \mathrm{1}+\mathrm{5}\mathit{\zeta }& \mathit{\zeta }\ge \mathrm{0}\end{array}\right.\end{array}$$

and

$$\begin{array}{}\text{(66)}& {\mathit{\varphi }}_{\mathrm{H}}\left(\mathit{\zeta }\right)=\left\{\begin{array}{ll}(\mathrm{1}-\mathrm{16}\mathit{\zeta }{)}^{-\mathrm{1}/\mathrm{2}}& \mathit{\zeta }<\mathrm{0}\\ \mathrm{1}+\mathrm{5}\mathit{\zeta }& \mathit{\zeta }\ge \mathrm{0},\end{array}\right.\end{array}$$

with ζ limited to the range −2 to 1. These functions integrate to give

$$\begin{array}{}\text{(67)}& \begin{array}{rl}& {\mathit{\psi }}_{\mathrm{m}}\left(\mathit{\zeta }\right)=\\ & \left\{\begin{array}{ll}\mathrm{2}\ln \left(\frac{\mathrm{1}+x}{\mathrm{2}}\right)+\ln \left(\frac{\mathrm{1}+x^{\mathrm{2}}}{\mathrm{2}}\right)-\mathrm{2}{\tan }^{-\mathrm{1}}x+\frac{\mathit{\pi }}{\mathrm{2}}& \mathit{\zeta }<\mathrm{0}\\ -\mathrm{5}\mathit{\zeta }& \mathit{\zeta }\ge \mathrm{0}\end{array}\right.\end{array}\end{array}$$

and

$$\begin{array}{}\text{(68)}& {\mathit{\psi }}_{\mathrm{H}}\left(\mathit{\zeta }\right)=\left\{\begin{array}{ll}\mathrm{2}\ln \left(\frac{\mathrm{1}+x^{\mathrm{2}}}{\mathrm{2}}\right)& \mathit{\zeta }<\mathrm{0}\\ -\mathrm{5}\mathit{\zeta }& \mathit{\zeta }\ge \mathrm{0},\end{array}\right.\end{array}$$

where $x=(\mathrm{1}-\mathrm{16}\mathit{\zeta }{)}^{\mathrm{1}/\mathrm{4}}$. The Obukhov length depends on fluxes, the fluxes depend on the stability functions and the stability functions depend on the Obukhov length, so stability adjustments are calculated iteratively starting from ζ=0.

2.5 Conducted heat fluxes

Snow and soil temperatures and liquid water fractions are simulated by a multi-layer heat conduction model (Essery, 2015). The conducted heat flux into the snow or ground surface is calculated as

$$\begin{array}{}\text{(69)}& G={\displaystyle \frac{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}{\mathrm{\Delta }{z}_{\mathrm{1}}}}(T_{\mathrm{s}}-T_{\mathrm{1}}),\end{array}$$

where λ₁ is the thermal conductivity of snow or soil and T₁ is the temperature of a surface layer of thickness Δz₁ (0.1 m by default).

2.6 Energy balance

2.6.1 Open areas

The surface energy balance,

$$\begin{array}{}\text{(70)}& f\left(T_{\mathrm{s}}\right)={\mathrm{LW}}_{\mathrm{s}}+{\mathrm{SW}}_{\mathrm{s}}-G-H-LE-L_{\mathrm{f}}M=\mathrm{0},\end{array}$$

with latent heat of fusion L_f and parametrizations for the fluxes, is a nonlinear equation for the unknown surface temperature and snowmelt rate M. The solution is first found with M=0. From an initial guess of temperature T_s0 and neglecting the complicated temperature dependence of r_a if stability adjustments are applied, a linear estimate of T_s is given by

$$\begin{array}{}\text{(71)}& \begin{array}{rl}{T}_{\mathrm{s}}& =T_{\mathrm{s}\mathrm{0}}-f\left(T_{\mathrm{s}\mathrm{0}}\right){\left({\displaystyle \frac{\mathrm{d}f}{\mathrm{d}{T}_{\mathrm{s}}}}\right)}^{-\mathrm{1}}=T_{\mathrm{s}\mathrm{0}}+f\left(T_{\mathrm{s}\mathrm{0}}\right)\\ & {\left[\mathrm{4}\mathit{\sigma }{T}_{\mathrm{s}\mathrm{0}}^{\mathrm{3}}+{\displaystyle \frac{\mathrm{2}{\mathit{\lambda }}_{\mathrm{1}}}{\mathrm{\Delta }{z}_{\mathrm{1}}}}+{\displaystyle \frac{\mathit{\rho }}{r_{\mathrm{a}}}}(c_p+L{D}_{\mathrm{s}}{\mathit{\psi }}_{\mathrm{s}})\right]}^{-\mathrm{1}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(72)}& D_{\mathrm{s}}={\left.{\displaystyle \frac{\mathrm{d}{q}_{\mathrm{sat}}}{\mathrm{d}T}}\right|}_{T=T_{\mathrm{s}\mathrm{0}}}={\displaystyle \frac{L{q}_{\mathrm{sat}}\left(T_{\mathrm{s}\mathrm{0}}\right)}{R_{\mathrm{wat}}{T}_{\mathrm{s}\mathrm{0}}^{\mathrm{2}}}}\end{array}$$

and R_wat is the gas constant for water vapour. A single evaluation of Eq. (71) gives an approximate solution, and repeated evaluations with T_s calculated in each iteration being used as T_s0 in the next is the Newton–Raphson method for solving Eq. (70). If this gives T_s>T_m and there is snow with ice mass I on the ground, iteration of Eq. (71) is repeated assuming that all of the snow melts and $M=I/\mathit{\delta }t$. If this gives T_s<T_m, the snow does not all melt and T_s=T_m is known; Eq. (70) is then solved instead for the unknown melt rate by substitution of T_s=T_m in the equations for the other fluxes.

2.6.2 Forest areas

For a one-layer canopy model, energy and water vapour mass conservation equations,

$$\begin{array}{}\text{(73)}& {\displaystyle }{f}_{\mathrm{1}}={\mathrm{LW}}_{\mathrm{s}}+{\mathrm{SW}}_{\mathrm{s}}-G-H_{\mathrm{s}}-{\mathrm{LE}}_{\mathrm{s}}-L_{\mathrm{f}}M=\mathrm{0},\text{(74)}& {\displaystyle }{f}_{\mathrm{2}}={\mathrm{LW}}_{\mathrm{v}}+{\mathrm{SW}}_{\mathrm{v}}-H_{\mathrm{v}}-{\mathrm{LE}}_{\mathrm{v}}-C_{\mathrm{v}}{\displaystyle \frac{\mathrm{d}{T}_{\mathrm{v}}}{\mathrm{d}t}}=\mathrm{0},\text{(75)}& {\displaystyle }{f}_{\mathrm{3}}=(H-H_{\mathrm{s}}-H_{\mathrm{v}})/\left(\mathit{\rho }{c}_p\right)=\mathrm{0},\end{array}$$

and

$$\begin{array}{}\text{(76)}& f_{\mathrm{4}}=(E-E_{\mathrm{s}}-E_{\mathrm{v}})/\mathit{\rho }=\mathrm{0},\end{array}$$

with parametrizations for the fluxes, form a set of four nonlinear equations with four unknowns: q_c, T_c, T_v, and either T_s or M. Writing vectors $\mathit{f}=(f_{\mathrm{1}},f_{\mathrm{2}},f_{\mathrm{3}},f_{\mathrm{4}}{)}^T$ and $\mathit{x}=(T_{\mathrm{s}},q_{\mathrm{c}},T_{\mathrm{c}},T_{\mathrm{v}}{)}^T$, a solution without melt is first found by iterating

$$\begin{array}{}\text{(77)}& \mathit{x}={\mathit{x}}_{\mathrm{0}}-{\mathbf{J}}^{-\mathrm{1}}\mathit{f}\left({\mathit{x}}_{\mathrm{0}}\right),\end{array}$$

where J is the Jacobian matrix of f with elements

$$\begin{array}{}\text{(78)}& J_{ij}={\displaystyle \frac{\partial f_i}{\partial x_j}}\end{array}$$

given in the full documentation distributed with the FSM2 code. Equation (77) is implemented by solving

$$\begin{array}{}\text{(79)}& \mathbf{J}(\mathit{x}-{\mathit{x}}_{\mathrm{0}})=-\mathit{f}\left({\mathit{x}}_{\mathrm{0}}\right)\end{array}$$

numerically and iterating to find x. If the solution has T_s>T_m and there is snow with ice mass I on the ground, Eq. (79) is iterated again with $M=I/\mathit{\delta }t$, assuming that all of the snow melts in the time step. If this gives T_s<T_m, the snow does not all melt; the surface temperature is then known and Eq. (79) is solved to find the unknown melt rate.

There are seven energy and water vapour mass conservation equations for the surface and two canopy layers and seven unknowns in a vector $\mathit{x}=(T_{\mathrm{s}},q_{\mathrm{c},\mathrm{1}},T_{\mathrm{c},\mathrm{1}},T_{\mathrm{v},\mathrm{1}},q_{\mathrm{c},\mathrm{2}},T_{\mathrm{c},\mathrm{2}},T_{\mathrm{v},\mathrm{2}}{)}^T$. The conservation equations and the elements of the 7×7 Jacobian matrix are, again, listed in the FSM2 documentation. Solutions are found in the same way as for the one-layer canopy.

2.7 Canopy snow

Early land surface models used the same interception capacities for liquid water and snow held on vegetation (e.g. the canopy capacity per unit VAI was 0.2 kg m⁻² in CLASS prior to version 3.1 and 0.1 kg m⁻² in CLM prior to version 5.0), but measured canopy snow loads can actually be much higher. Subsequently, many models have adopted the representation of snow interception and unloading developed from observations by Hedstrom and Pomeroy (1998). Lundquist et al. (2021) noted that snow interception in global models is still based on a few geographically limited observations.

In FSM2, a forest canopy intercepts a fraction of falling snow up to a maximum S_c=S_ΛΛ, with S_Λ=4.4 kg m⁻² by default (Essery et al., 2003). The fraction of the canopy covered with snow, which is required for canopy albedo and sublimation calculations, is parametrized as

$$\begin{array}{}\text{(80)}& f_{\mathrm{cs}}={\left({\displaystyle \frac{S_{\mathrm{v}}}{S_{\mathrm{c}}}}\right)}^{\mathrm{2}/\mathrm{3}}\end{array}$$

for a canopy layer with intercepted snow mass S_v. The $\mathrm{2}/\mathrm{3}$ exponent here is quoted in many papers without citation or explanation. In fact, it was introduced by Deardorff (1978), who proposed it as a compromise between values of 0 and 1 that would make evaporation of dew from vegetation too fast and too slow, respectively.

The mass balance equation for snow in a canopy layer with interception rate I_v, sublimation rate E_v, melt rate M_v, and unloading rate U_v is

$$\begin{array}{}\text{(81)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{v}}}{\mathrm{d}t}}=I_{\mathrm{v}}-E_{\mathrm{v}}-M_{\mathrm{v}}-U_{\mathrm{v}}.\end{array}$$

If vegetation temperature T_v0 calculated without melt exceeds T_m while intercepted snow remains, the canopy snowmelt rate is

$$\begin{array}{}\text{(82)}& M_{\mathrm{v}}=min\left[{\displaystyle \frac{C_{\mathrm{v}}}{L_{\mathrm{f}}\mathit{\delta }t}}(T_{\mathrm{v}\mathrm{0}}-T_{\mathrm{m}}),{\displaystyle \frac{S_{\mathrm{v}}}{\mathit{\delta }t}}\right]\end{array}$$

and the vegetation temperature is reset to

$$\begin{array}{}\text{(83)}& T_{\mathrm{v}}=T_{\mathrm{v}\mathrm{0}}-{\displaystyle \frac{L_{\mathrm{f}}{M}_{\mathrm{v}}\mathit{\delta }t}{C_{\mathrm{v}}}}.\end{array}$$

Snowfall that is not intercepted is added to the snow on the ground with the same fresh snow density as for snow falling in open areas, and snow unloading from the canopy is added with the same density as snow already on the ground; this is a crude approximation, but model development is limited by a lack of measurements of the evolution of intercepted snow properties (Bouchard et al., 2024). Meltwater dripping from the canopy and rain are added to the snow as liquid water. Interception, unloading, and melt are calculated for both layers in the two-layer canopy model; the lower layer can intercept throughfall of snow but not unloading or drip from the upper layer.

2.7.1 Canopy interception options

CLM and MATSIRO intercept constant fractions of snowfall until the canopy capacity is reached. This linear option is implemented in FSM2 as

$$\begin{array}{}\text{(84)}& I_{\mathrm{v}}=min\left[f_{\mathrm{v}}{S}_{\mathrm{f}},{\displaystyle \frac{(S_{\mathrm{c}}-S_{\mathrm{v}})}{\mathit{\delta }t}}\right]\end{array}$$

for snow falling at the rate S_f. CLASS, ISBA, JULES, and VISA use the Hedstrom and Pomeroy (1998) interception rate model

$$\begin{array}{}\text{(85)}& I_{\mathrm{v}}={\displaystyle \frac{(S_{\mathrm{c}}-S_{\mathrm{v}})}{\mathit{\delta }t}}\left[\mathrm{1}-\exp \left(-{\displaystyle \frac{f_{\mathrm{v}}{S}_{\mathrm{f}}\mathit{\delta }t}{S_{\mathrm{c}}}}\right)\right],\end{array}$$

which gives the same interception rate as Eq. (84) for an initially snow-free canopy but approaches S_c more slowly. Equations (84) or (85) are applied in both layers of the two-layer canopy model, with interception in the upper layer subtracted from snowfall reaching the lower layer. Because of the non-linearity in Eq. (85), one-layer and two-layer representations of the same canopy density have different interception efficiencies.

2.7.2 Canopy unloading options

Canopy snow loads decrease exponentially with time after snowfall in CLASS, JULES, and MEB. This time-dependent option is implemented in FSM2 along with increased unloading when canopy snow is melting as

$$\begin{array}{}\text{(86)}& U_{\mathrm{v}}={\displaystyle \frac{S_{\mathrm{v}}}{{\mathit{\tau }}_{\mathrm{u}}}}+m_{\mathrm{u}}{M}_{\mathrm{v}},\end{array}$$

with τ_u=10 d (Bartlett and Verseghy, 2015) and m_u=0.4 by default (Storck et al., 2002).

CLM5 and VISA have unloading rates that depend on canopy temperature and wind speed. This temperature- and wind-dependent unloading option is implemented in FSM2 as

$$\begin{array}{}\text{(87)}& U_{\mathrm{v}}=\left[{\displaystyle \frac{\mathrm{1}}{c_T}}max(T_{\mathrm{v}}-\mathrm{270.15},\mathrm{0})+{\displaystyle \frac{U_{\mathrm{a}}}{c_U}}\right]S_{\mathrm{v}},\end{array}$$

with $c_T=\mathrm{1.87}\times {\mathrm{10}}^{\mathrm{5}}$ K s and $c_U=\mathrm{1.56}\times {\mathrm{10}}^{\mathrm{5}}$ m by default (Roesch et al., 2001). With the default parameter values, unloading rates from Eq. (87) exceed rates from Eq. (86) whenever the canopy temperature is above −2.8 °C or the wind speed is above 0.2 m s⁻¹.

3 Model test results

The following comparisons of FSM2 results with observations are presented as demonstrations of model behaviour rather than rigorous evaluations of the model. Simulations of sub-canopy shortwave radiation, longwave radiation, and wind speeds are first compared with observations at sites in Finland, Sweden, and Switzerland. A complete simulation of snow mass and energy balance is then compared with a year of observations at a site in Switzerland. Finally, FSM2 simulations are compared with observations and Land-Surface, Snow and Soil moisture Model Intercomparison Project (LS3MIP) simulations of Northern Hemisphere snow cover extent, snow mass, and albedo.

3.1 Sub-canopy radiation

Radiation was measured with arrays of 10 shortwave radiometers and four longwave radiometers under sparse forest canopies over periods of between 4 and 23 d during March and April 2011 at Abisko, Sweden (68.3° N, 18.8° E) and during March and April 2012 at Sodankylä, Finland (67.4° N, 26.6° E). Measurements were made in five stands described by Reid et al. (2013) at each site: leafless birch stands with average sky view fractions determined from hemispherical photography between 0.59 and 0.9 at Abisko and pine, spruce, and mixed stands with average sky view fractions between 0.41 and 0.72 at Sodankylä. There was snow on the ground, but the canopy was snow-free during all of the measurement periods. For consistency in the model, vegetation area indices were obtained from sky view fractions by inverting Eq. (14). Above-canopy meteorological driving datasets were constructed using measurements from the Abisko Scientific Research Station and the Finnish Meteorological Institute Arctic Space Centre at Sodankylä.

Figure 5 shows average sub-canopy shortwave radiation and shortwave radiation transmission calculated as average sub-canopy radiation divided by average above-canopy radiation over the measurement periods. Simulations are shown with Beer's law and two-stream radiative transfer options; the matrix formulation of the radiative transfer equations ensures that one- and two-layer canopy models give the same results if the albedos of the layers are the same. As was seen in Fig. 3, diffuse transmission is slightly higher and direct-beam transmission is lower for snow-free canopies if Beer's law is used. Averaged over the measurement periods for each of the stands, the simulated transmission in Fig. 5b is lower for Beer's law, particularly for measurement periods during which lower fractions of the cumulated incoming shortwave radiation were diffuse (these fractions were between 34 % and 59 % over the measurement periods).

Figure 5 Beer's law (open circles) and two-stream approximation (closed circles) simulations of (a) average sub-canopy downward shortwave radiation and (b) canopy shortwave transmission, compared with measurements in 10 forest stands. Simulations with one and two canopy layers are indistinguishable.

Figure 6 Simulations with one canopy layer (open circles) and two canopy layers (closed circles) of (a) average sub-canopy downward longwave radiation and (b) canopy longwave enhancement, compared with measurements in 10 forest stands.

In analogy with shortwave radiation transmission, dividing average sub-canopy longwave radiation by average above-canopy longwave radiation gives a canopy longwave enhancement factor. This is generally greater than one because canopy elements have higher emissivities than the atmosphere, particularly when the sky is clear (Rutter et al., 2023). The canopy longwave enhancement emphasizes differences in model canopy temperature rather than differences in above-canopy longwave radiation for the different measurement periods in Fig. 6. The upper-canopy layer shades the lower layer from shortwave radiation by day and traps longwave radiation by night in a two-layer canopy model (Gouttevin et al., 2015), but average differences between FSM2 simulations with one and two layers are small for these sparse canopies.

3.2 Sub-canopy wind speed

Wind speed was measured by single 2D sonic anemometers at 2 m height in open areas and under forest canopies in four pine stands (8 February to 5 March 2018) and 14 spruce stands (10 January to 25 April 2018 and 24 January to 13 April 2019) near Davos, Switzerland (46.8° N, 9.9° E) and four pine stands near Sodankylä (17 to 29 April 2019), with Λ ranging from 1 to 3.4 (Mazzotti et al., 2020a). Figure 7 compares simulated forest wind speeds and ratios of forest to open wind speeds averaged over the measurement periods, which varied from 3 to 58 d in length. The observations may be influenced by forest edge effects, and there is a large degree of scatter in the simulations, but the simulated ratios have a moderate correlation (0.64) with the observations and a moderate root mean square error (0.1). Equation (58) gives a continuous vertical profile of wind speeds without discretization, so results for one- and two-layer canopy models do not differ.

Figure 7 Simulations of (a) average sub-canopy wind speeds and (b) ratios between sub-canopy and open wind speeds, compared with measurements in 22 forest stands. Simulations with one and two canopy layers are indistinguishable.

3.3 Snow simulations at a site

A site has to be well characterized, well instrumented, and well maintained to provide direct measurements of all of the inputs required by energy balance models. Several such sites have been used in the Snow Model Intercomparison Project (SnowMIP; Etchevers et al., 2004; Essery et al., 2009; Menard et al., 2019) for evaluation of snow models. Forest and open meadow sites at Alptal, Switzerland (Stähli and Gustafsson, 2006; Stähli et al., 2009), that were used in SnowMIP2 were also used by Gouttevin et al. (2015) and are used again here to demonstrate the performance of FSM2. The forest stand is dominated by spruce and fir with typical heights of 25 m and an average VAI of 3.96. The mild climate can allow snow cover to appear and disappear several times over the winter; such conditions are known to be challenging for snow modelling (Essery et al., 2009). FSM2 was run for the winter of 2004–2005 at Alptal, with all driving data except precipitation measured above the forest canopy on a 35 m high mast. Precipitation was measured with a sheltered gauge in the meadow and gauge-corrected by scaling snowfall to match snow accumulation in the meadow between snow-free conditions observed on 13 December 2004 and the maximum snow mass of 352 kg m⁻² observed on 14 March 2005. Model parameters α_c0, α_Λ0, and α₀ were adjusted to match measured snow-free albedos above the forest and the meadow, but all other parameters were left at default values to focus on differences in simulations due to differences in model options. Figure 8 compares simulations with observations from Stähli et al. (2009) of average snow mass measured weekly along 30 m transects in the forest and in the meadow, albedo and outgoing (upward) longwave radiation measured at fixed points above the forest and in the meadow, shortwave transmission and incoming (downward) longwave radiation measured with radiometers moving along a 10 m horizontal rail under the canopy, and a subjective estimate of canopy snow load made by an observer from 0 (snow-free canopy) to 8 (maximum possible snow interception), scaled to the model's canopy capacity. The albedo, transmission, and incoming longwave radiation measurements have been filtered to remove periods when it is suspected that there was snow on the upwards-facing radiometers. Lower snow mass, lower albedo, and higher outgoing longwave radiation for the forest than for the meadow are apparent in both observations and simulations.

Figure 8 Simulations (lines) and observations (points) for the Alptal forest and meadow sites in 2004–2005. Albedo, transmission, and longwave radiation fluxes are daily values. Subjective canopy snow load observations are scaled to the model's canopy capacity (dashed line in panel b).

The same model options were used for snow on the ground in the meadow and the forest, corresponding to FSM configuration 31: prognostic snow albedo, variable thermal conductivity, prognostic snow density, stability adjustment of the turbulent exchange coefficient, and prognostic liquid water content. Because cumulated snowfall in the model driving data was forced to match the observed maximum snow mass in the meadow, it is not surprising that the simulation matches the snow accumulation well, but the model also matches observed snowmelt well, including short periods with shallow snow cover in November 2004 and April 2005. Only a shallow snow cover is required to increase the meadow albedo from below 0.2 to above 0.7.

There are 16 simulations for the forest: with linear or nonlinear interception efficiency, with one or two canopy layers, with Beer's law or two-stream canopy radiative transfer, and with time- and melt-dependent or temperature- and wind-dependent unloading. The forest simulations have a 57 kg m⁻² range in maximum snow mass and a 46 d range in duration of snow cover (Fig. 9 and Table 1). Although large compared with plausible observation errors, these ranges are smaller than the spread seen in comparisons of forest snow simulations by different models (Rutter et al., 2009, and Appendix A); Essery et al. (2009) suggested that this spread was largely due to uncertain parameter selections for highly simplified canopy models. The forest canopy is dense enough that snow on the ground has very little influence on above-canopy albedo (Fig. 8c). The albedo increases when there is intercepted snow in the canopy but remains below 0.3. Agreement between the durations of periods with elevated albedo in observations and simulations suggests that the simulated persistence of snow in the canopy is realistic. Simulated transmission of shortwave radiation through the canopy (Fig. 8e) increases when there is intercepted snow in the canopy but is lower than observed; the transmission could be increased with little impact on simulated snow masses by decreasing the parameter k_ext. In the lower layer of a two-layer canopy model, daytime heating by shortwave radiation and nighttime cooling by longwave radiation are reduced by the shelter of the upper layer (Gouttevin et al., 2015; Todt et al., 2018). This reduces the diurnal range in sub-canopy longwave radiation, but differences in simulations of daily-average longwave radiation (Fig. 8d and f) are small.

There are too many lines in Fig. 8a to identify the individual canopy model configurations, so Table 1 and Fig. 9 give the maximum sub-canopy snow mass, the duration of snow cover on the ground, and the fraction of total snowfall sublimating for each simulation. Increasing sublimation is strongly related to decreasing snow mass (correlation $r=-\mathrm{0.93}$), and decreasing snow mass is strongly related to decreasing snow cover duration (r=0.96). The unloading option has the largest influence on snow on the ground; snow unloads from the canopy faster with the temperature- and wind-dependent unloading option and accumulates on the ground, where it is sheltered from wind and sublimation. Less snowfall is intercepted by the nonlinear interception option as the canopy load increases, so this option also increases the mass of snow on the ground. Differences in transmission of shortwave radiation through the canopy between the Beer's law and two-stream radiative transfer options depend on canopy snow and sky conditions (Fig. 3b and d), but these differences have little influence on snow beneath the canopy because the transmission is always low for the dense Alptal forest canopy. The number of model canopy layers has complex influences on snow simulations. The upper layer in a two-layer model intercepts more snowfall than the lower layer, and that snow is exposed to higher wind speeds and shortwave radiation. Snow in the lower layer is sheltered, but the overall effect is for slightly more snow sublimation. The shading of the lower layer, however, decreases daytime sub-canopy longwave radiation, which reduces mid-winter melt and delays the final disappearance of snow on the ground in spring.

Table 1 Maximum sub-canopy snow mass, duration of snow cover on the ground, and fraction of total snowfall sublimating in the 16 forest simulations in Fig. 8a with every possible combination of linear or nonlinear snowfall interception, one or two canopy layers, Beer's law or two-stream canopy radiative transfer, and time- and melt-dependent (t−M) or temperature- and wind-dependent (T−U) canopy snow unloading.

Figure 9 Maximum sub-canopy snow mass, duration of snow cover on the ground, and fraction of total snowfall sublimating in the 16 Alptal forest simulations numbered in Table 1.

3.3.1 Sensitivity to canopy density

The amounts of snowfall reaching the ground in open and forest sites differ because some of the intercepted snow in the canopy sublimates, and meltwater dripping from snow in the canopy drains from the snow on the ground if it does not refreeze. Melt rates on the ground differ between sites because the canopy modifies the shortwave radiation, longwave radiation, and turbulent heat fluxes in the surface energy balance. All of these differences are influenced by the canopy density, which is represented by the vegetation area index in FSM2. Figure 10 explores variations with VAI in simulations using the Alptal 2004–2005 meteorology.

Figure 10a shows fractions of snowfall that sublimate from the ground and canopy combined. Light winds (below 4 m s⁻¹ for 97 % of hours) give low sublimation for the open meadow site (Λ=0), but sublimation from the canopy increases as the density increases. By design in FSM2, simulations are continuous as Λ→0 for sparse canopies and are independent of canopy density for dense canopies with Λ≳4. Differences between simulations increase as canopy density increases and level out for high canopy densities, limited by energy availability. Reference to Table 1 again allows identification of the individual model configurations in Fig. 10a; simulations with the time- and melt-dependent unloading option consistently have the highest sublimation because snow remains exposed in the canopy for longer.

Figure 10 Sensitivity to canopy density in simulations with the 16 canopy model configurations and the Alptal 2004–2005 meteorology. (a) Fractions of total snowfall sublimating. (b) Contributions of net shortwave radiation, net longwave radiation, and sensible heat fluxes to energy for melting snow on the ground.

The contributions of energy by net shortwave radiation, net longwave radiation, and sensible heat fluxes to melt snow on the ground in simulations with varying canopy density are shown in Fig. 10b. All 16 canopy model configurations give similar trends with canopy density. Simulated melt at the Alptal meadow site is dominated by shortwave radiation, with a small contribution from sensible heat fluxes, and net longwave radiation is a small loss of energy for snowmelt. As the canopy density increases and sky view under the canopy decreases, the net shortwave radiation decreases and the net longwave radiation increases to become the dominant source for melt energy under dense canopies. Sensible heat fluxes first increase as canopy air space temperatures increase because of canopy heating and then decrease as sub-canopy wind speeds decrease with increasing canopy density. Latent and ground heat fluxes (not shown) each contribute less than 10 % of the melt energy beneath canopies.

3.4 Northern Hemisphere snow simulations

Simulation of snow on a grid covering an area requires coupling with an atmospheric model or distributed driving data that are not directly available from measurements. Meteorological reanalyses can be used for simulations over large areas at coarse resolutions; for example, Brun et al. (2013) used ERA-Interim reanalyses to drive the Crocus snow model over northern Eurasia. Here, the performance of FSM2 for simulating Northern Hemisphere seasonal snow cover at 0.5° resolution is demonstrated with 2000–2010 driving data from the GSWP3 bias-corrected reanalysis (Kim, 2017), which was previously used in LS3MIP (van den Hurk et al., 2016). Canopy heights, vegetation area indices, forest fractions, and snow-free albedos for FSM2 are taken from global maps developed by Lawrence and Chase (2007) for CLM; deciduous and evergreen forest fractions are shown in Fig. 11. Separate FSM2 simulations with parameters for evergreen forest, deciduous forest, and unforested land are combined according to their fractions to give the averaged seasonal cycles of Northern Hemisphere snow area and mass shown in Fig. 12 for the 16 canopy model configurations listed in Table 1. The same model configuration as in Sect. 3.3 was used for unforested land, and snow interception was set to zero for leafless deciduous forest canopies, so differences between simulations are dominated by areas with evergreen forests.

Figure 12 compares FSM2 simulations with the nine models that submitted results for LS3MIP and estimates from multi-dataset historical snow extent and snow mass time series (Mudryk et al., 2020). One of the LS3MIP models follows the estimated snow mass closely but underestimates snow extent. FSM2 and the other LS3MIP models give very similar snow masses while persistent snow cover accumulates at high latitudes through October and November. Thereafter, the models retain different amounts of ephemeral snow mass and spread out. FSM2 tends to overestimate snow mass and underestimate the peak snow area in comparison with the historical time series but lies within the ranges of the LS3MIP models; Mudryk et al. (2020) noted that the observational estimates of hemispheric snow mass are likely to be biased low because of underestimation in mountain regions. Despite noticeable local differences in snow mass on the ground under forests demonstrated in Sect. 3.3, the spread in hemispheric averages for FSM2 simulations in Fig. 12 is small. This is because evergreen needleleaf forests only cover 8 % of the Northern Hemisphere land area in the dataset of Lawrence and Chase (2007) and around 20 % of the area with seasonal snow cover in FSM2 simulations.

Figure 11 Northern Hemisphere deciduous and evergreen forest fractions from Lawrence and Chase (2007).

Figure 12 Northern Hemisphere seasonal snow area (a) and mass (b) averaged over 2000–2010 from 16 FSM2 simulations (black lines), nine LS3MIP models (blue lines), and estimates from multi-dataset historical time series (circles).

As an example of the spatial distribution of Northern Hemisphere snow cover, Fig. 13 compares FSM2 with binary snow cover information at 4 km resolution from the Interactive Multi-sensor Snow and Ice Mapping Service (IMS; US National Ice Center, 2008) and white sky albedo at 0.05° resolution from the MODIS MCD43C3 dataset (Schaaf and Wang, 2021) aggregated to 0.5° resolution on 1 March 2010. The good match in hemispheric snow cover extent at this time of year has some compensation between underestimates and overestimates at the southern limits of snow cover (Fig. A3 shows similar snow cover difference maps for the LS3MIP models). Albedo differences are largest where FSM2 has errors in snow cover fraction. FSM2 and MODIS both show dark bands across the continents where the albedo of snow is masked by boreal forests, but the FSM2 albedo is generally higher. Evaluation and optimization of FSM2 for hemispheric snow simulations would clearly require much larger samples of observations and simulations, but these preliminary results are encouraging.

Figure 13 Northern Hemisphere snow cover fraction and albedo maps for 1 March 2010 from FSM2, IMS, and MODIS MCD43C3. The FSM2 simulation has two canopy layers, two-stream radiative transfer, nonlinear canopy snow interception, and temperature- and wind-dependent unloading.

4 Discussion

Despite long-standing interest, uncertainty remains in how best to represent canopy and sub-canopy snow processes in models. In contrast with FSM2, Gouttevin et al. (2015) made a number of different design decisions when implementing a canopy model in SNOWPACK. For a two-layer model, SNOWPACK conceptualizes the canopy as having an upper leaf layer and a lower trunk layer, without multiple reflections of shortwave radiation in the upper layer or interception of snow in the lower layer. Considering that turbulent transport within vegetation canopies is still poorly understood, Gouttevin et al. (2015) followed Blyth et al. (1999) in using logarithmic wind profiles even within canopies and making empirical adjustments to surface aerodynamic resistances beneath canopies for simplicity. The two vegetation layers in SNOWPACK are coupled to a single canopy air space ($T_{\mathrm{c},\mathrm{1}}=T_{\mathrm{c},\mathrm{2}}$, $r_{\mathrm{v},\mathrm{1}}=r_{\mathrm{v},\mathrm{2}}$ and r_c=0 in Fig. 2), and a minimum heat exchange coefficient for windless conditions is included in the parametrizations of fluxes between the canopy air space and the atmosphere. The aerodynamic resistance between the canopy and the canopy air space in SNOWPACK (Eq. 32 in Gouttevin et al., 2015) does not scale with leaf area and can be more than 2 orders of magnitude smaller than the same resistance in FSM2 (Eq. 63). Using two-layer canopy options decreases the diurnal ranges of sub-canopy downward longwave radiation in both FSM2 and SNOWPACK. In SNOWPACK, Gouttevin et al. (2015) show in their Fig. 3 that this reduction is dominated by an increase in nighttime minima, and daily averages are increased. Heating of the lower layer on clear days and cooling on clear nights are both reduced in FSM2, and differences in daily averages are small (Fig. 8 here). Differences in simulations of maximum snow mass on the ground at Alptal in March 2005 between one-layer and two-layer canopy models are smaller in FSM2 (up to 16 kg m⁻² in Table 1) than in SNOWPACK (54 kg m⁻² in Fig. 4 of Gouttevin et al., 2015).

Unloading of canopy snow does not occur in SNOWPACK until the intercepted snow load exceeds the canopy capacity (Eq. 5 in Gouttevin et al., 2015), whereas unloading is continuous in FSM2 (Eqs. 86 or 87). The canopy snow unloading parametrizations reviewed by Lundquist et al. (2021) gave the largest differences in simulations of snow on the ground at Alptal when coupled in the full mass and energy balances of FSM2. Lundquist et al. (2021) found little experimental evidence for a maximum snow interception capacity in published datasets. In fact, there is also little evidence for a maximum capacity in the Alptal forest simulations because the model's capacity is rarely reached (Fig. 8b). Rain falling on snow is a rare phenomenon occurring most frequently in the Arctic and northern maritime regions that can have large impacts (Cohen et al., 2015); rain does not interact with canopy snow in FSM2, but this has been investigated in SNOWPACK by Bouchard et al. (2024).

The two-layer and two-stream models in FSM2 are the more physical canopy structure and radiative transfer options, but they have relatively little influence on simulations of sub-canopy snow for the Alptal test case. The empirical canopy snow interception and unloading options make larger differences but have not been widely tested in different climates and forest types. FSM2 is formulated to limit differences between one-layer and two-layer canopy models for reasons that are purely related to the methods chosen for numerical solution of the mass and energy balance equations. If the albedos and temperatures of both layers in a two-layer model are equal, the sub-canopy shortwave radiation and longwave radiation will be the same as in a one-layer model. The Hedstrom and Pomeroy (1998) snow interception, however, was developed as a bulk canopy model. The nonlinearity of Eq. (85) means that a two-layer canopy model has a slightly higher interception efficiency than a one-layer model with the same canopy capacity. The vertical distribution of snow in the canopy is important for sublimation because snow higher in the canopy is exposed to higher wind speeds and higher solar radiation. Bonan et al. (2021) argued that 5 to 10 canopy layers are necessary to model turbulent and radiative fluxes but did not consider interception. Terrestrial laser scanning can now be used to make measurements of canopy loading (Russell et al., 2021) that might be useful in developing an explicitly multilayer interception model.

Qu and Hall (2014) reported a large spread between climate models for the albedo of snow-covered land in boreal forest regions, with implications for simulations of snow–albedo feedback. Different canopy radiative transfer parametrizations, however, need not result in large differences in masking of snow albedo by forests. FSM2 parameters have been adjusted to give similar canopy albedos from Beer's law and the two-stream approximation; transmission of shortwave radiation through the canopy differs, but this has little influence on snowmelt beneath dense canopies. Variations in canopy parameters, VAI maps, and alternative canopy structure metrics have not been explored here but can give large variations in canopy albedo (Essery, 2013; Malle et al., 2021).

5 Outlook

The addition of canopy model options enables the use of FSM2 for investigating snow–forest interactions and for large-scale simulations including forested regions. During the course of its development, FSM2 has been coupled with the Snow Microwave Radiative Transfer model (SMRT; Sandells et al., 2017), the Multiple Snow Data Assimilation System (MuSA; Alonso-González et al., 2022), and the Swiss operational snow-hydrological model system (OSHD; Mott et al., 2023). It has been adapted for metre-resolution simulation of snow in discontinuous forests (Mazzotti et al., 2020a, b, 2023) by removing the assumption of homogeneous canopy cover for radiative transfer calculations. Implementations of FSM2 in OSHD (Quéno et al., 2024) and the Canadian Hydrological Model (CHM; Marsh et al., 2020) add representations of horizontal snow redistribution at snowdrift-resolving scales. Parametrizations for interception of snow by deciduous conifers such as larch and trapping of drifting snow by tundra shrubs would be other useful additions. FSM2 is now being coupled with the Open Global Glacier Model (OGGM; Maussion et al., 2019) for physically based mass balance modelling of glaciers; representations of firn compaction (Lundin et al., 2017) and debris cover (Reid et al., 2012) will be added for this. A limitation of FSM2 is that it does not have a dynamic representation of soil moisture; it has to be coupled with groundwater and routing models for many hydrological applications.

In a comparison with observations on small scales at sites with varying climatic conditions, FSM2 captures the differences in snow mass, snow cover duration, and albedo between open and forested sites. For Northern Hemisphere seasonal snow area and mass simulations, FSM2 lies within the range of land surface schemes from state-of-the-art Earth system models. Because FSM2 is intended for snow research and does not attempt to include all the land processes required in modern Earth system models, it is more compact and easier to use than these land surface schemes. There are 3186 lines of code in the FSM2 source directory, compared with more than 190 000 for CLM5.0, and the FSM2 code compiles in under 5 s. The Northern Hemisphere simulations presented here took about 7 min of CPU time to run 50 411 land points for 2920 time steps per year on a single 1.5 GHz processor. No communication is required between points, so the code is trivially parallelizable. FSM2 offers opportunities for investigations requiring snow simulations on large grids, in large ensembles, or for long periods.

Appendix A: Snow model intercomparisons

Figure A1 Snow mass simulated with the 16 canopy configurations of FSM2 (black lines), simulated by the 33 models that participated in SnowMIP2 (grey lines) and measured at the four SnowMIP2 forest sites (red points). Snow mass measurements for the first winter at each site were shared with the SnowMIP2 participants to allow model calibration, but the models still produced a wide range of results (Rutter et al., 2009).

Figure A2 Snow mass simulated with the 16 canopy configurations of FSM2 (black lines), simulated by the 23 models that participated in ESM-SnowMIP (grey lines) and measured at the three ESM-SnowMIP forest sites (red points) (Menard et al., 2019). No data were provided for model calibration. The large underestimates in simulated snow mass at the aspen site in 2007–2008 resulted from erroneous input snowfall data.

Figure A3 Differences between simulations of Northern Hemisphere snow cover by the nine models that participated in LS3MIP (van den Hurk et al., 2016) and IMS on 1 March 2010.