Enhancing simulations of snowpack properties in land surface models with the Soil, Vegetation and Snow scheme v2.0 (SVS2)

Snow microstructure – characterized by density, grain size, grain shape and arrangement – fundamentally determines snowpack macroscopic properties. Despite this critical role, many land surface models (LSMs) lack explicit representation of snow microstructure. This limitation has become increasingly critical as future spaceborne missions for snow water equivalent measurement demand advanced modelling systems capable of accurately estimating snowpack properties, including microstructure, across diverse climatic and vegetation regions. The Soil Vegetation and Snow (SVS) LSM, used by Environment and Climate Change Canada for operational land surface and hydrological predictions, has been substantially upgraded to address these challenges. SVS version 2.0 (SVS2) incorporates the detailed multilayer Crocus snowpack model, enabling distinct simulations of snowpack evolution in both open terrain and forested areas within each grid cell. Crocus within SVS2 has been upgraded from its original alpine design with three major enhancements to handle Canada's varied snowpack conditions: an advanced albedo parameterization that accounts for spatial variability in light-absorbing particle deposition, new physical parameterizations tailored to Arctic snow characteristics, and a refined canopy model for forest environments. Significant improvements in simulations of near-surface density predictions are evident along a latitudinal transect from southern Quebec to the Canadian Arctic, while challenges remain in simulation of density and specific surface area in basal snow layers. SVS2 achieved substantial gains in snow melt-out timing accuracy, reducing prediction errors by over 50 % compared to the alpine Crocus version and surpassing two established snow reference products (ERA5-Land and ERA5-Crocus). These enhancements position SVS2 as a substantial advancement for future operational snow modeling applications across Canada.

1 Introduction

Snow cover plays a crucial role in Earth's climate system, influencing surface energy balance, water resources, ecosystem dynamics and ground thermal regime across vast regions of the planet (e.g., Armstrong and Brun, 2008; Flanner et al., 2011; Trujillo et al., 2012; Park et al., 2015; Sturm et al., 2017; Immerzeel et al., 2020). The role of snow is explained by unique physical properties such as its high albedo, the fraction of incident solar radiation reflected by the snow, and its low thermal conductivity governing heat transfer through the snowpack (Warren, 1982; Sturm et al., 1997; Armstrong and Brun, 2008). Snow physical properties are governed by the geometrical arrangement of individual components (air, ice, and possibly water and impurities) at the micrometer scale, referred to as snow microstructure (Fierz et al., 2009). For example, snow albedo in the near-infrared range strongly depends on snow specific surface area and the shape of ice particles (Domine et al., 2006; Robledano et al., 2023), while snow thermal conductivity depends on snow density and anisotropy (Yen, 1981; Calonne et al., 2011; Leinss et al., 2020). Snow microstructure also controls liquid water content in the snowpack (Yamaguchi et al., 2010) and its microwave scattering properties (Picard et al., 2022). Snow microstructural properties can strongly vary with depth, and control penetration of solar radiation, conduction of heat, and percolation of meltwater through the snowpack (e.g. Wever et al., 2016; Domine et al., 2016b; Barrere et al., 2017).

Land surface models (LSMs) serve as the terrestrial component of Earth System Models, Numerical Weather Prediction (NWP) systems, and hydrological models, representing complex interactions between the atmosphere and land surface processes (e.g., Overgaard et al., 2006; Fisher and Koven, 2020; Blyth et al., 2021). Within all of these models, snowpack schemes are essential components that simulate accumulation, evolution, and melt of seasonal snow cover (e.g., Boone and Etchevers, 2001; Niu et al., 2011; Vionnet et al., 2012; Arduini et al., 2019). They solve the energy and mass balance equations using approaches of varying complexity (Essery et al., 2013; Lafaysse et al., 2017). The representation of snowpack processes in land surface models presents a fundamental trade-off between computational efficiency and physical realism. One or two layer snowpack schemes represent a simplified snowpack and only simulate its bulk properties. Such simplified models are used operationally in certain climate models (Melton et al., 2020), NWP systems (Douville et al., 1995; Bélair et al., 2003), and hydrological models (Marks et al., 1999; Leonardini et al., 2021; Pomeroy et al., 2022).

Multi-layer snowpack schemes of intermediate complexity have been developed within LSMs to improve the representation of the snowpack internal properties (e.g., Boone and Etchevers, 2001; Niu et al., 2011; Essery et al., 2013; Decharme et al., 2016; Arduini et al., 2019; Lawrence et al., 2019; Pflug et al., 2019). These schemes vertically discretize the snowpack with a prescribed number of layers (up to 12) of fixed maximum thickness and account for some internal processes, such as snow settling, water percolation and refreezing. Cristea et al. (2022) have shown that such multi-layer schemes improve the representation of heat transport and mass storage within the snowpack. By incorporating multi-layer snowpack representations into their respective models, Saha et al. (2017) and Arduini et al. (2019) demonstrated enhanced accuracy in simulating both 2 m air temperature and snow depth across snow-covered areas. These improvements resulted from better representation of surface snow thermal conductivity and albedo, as well as more realistic snowpack thermal regime. Schemes of intermediate complexity do not generally include an explicit representation of snow microstructure and many snow properties are derived from snow density (Boone and Etchevers, 2001; Essery et al., 2013).

Among the 21 snowpack models used in the Earth System Model – Snow Model Intercomparison Project (ESM-Snow MIP, Menard et al., 2021), only three of them (Crocus (Brun et al., 1992; Vionnet et al., 2012), SNOWPACK (Bartelt and Lehning, 2002), and the Snow Metamorphism and Albedo Process model (SMAP Niwano et al., 2012)) explicitly simulate the evolution of snow microstructure, through snow metamorphism. These detailed multi-layer snowpack schemes simulate evolution of vertical layering without pre-defined assumptions about this layering for a given total snow height. These schemes can be used to simulate formation of distinct layers, such as surface hoar (Horton and Jamieson, 2016) or ice layers resulting from freezing rain (Quéno et al., 2018). Despite their ability to simulate advanced snow processes, detailed snowpack schemes are rarely implemented in LSMs nor used over large extents due to concerns about computational costs and associated data storage. For example, SNOWPACK is not implemented in any LSM with the exception of the CRYOWRF coupled snow/atmosphere model (Sharma et al., 2023). On the other hand, Crocus is implemented as a snow option in the SURFace EXternalisée (SURFEX) modelling platform (Vionnet et al., 2012; Masson et al., 2013), but its operational use is restricted to avalanche hazard forecasting over relatively small regions (Lafaysse et al., 2013). Large scale applications of Crocus within SURFEX (Brun et al., 2013; Ramos Buarque et al., 2025) or within the Land Information System (LIS, Navari et al., 2024) have nonetheless shown potential to provide accurate estimation of bulk snowpack properties (depth, density, snow water equivalent, SWE) across large spatial areas.

Figure 1 Land surface tiling used in SVS2 and overview of main snow-related components

Development of future spaceborne missions to retrieve SWE (Derksen et al., 2021; Tsang et al., 2021) requires state-of-the-art modelling systems to provide accurate estimation of snowpack properties (including microstructure) across large areas with contrasting climate and vegetation. Snow microstructural properties are required to inform SWE retrieval (King et al., 2018; Rutter et al., 2019; Singh et al., 2024) and data assimilation of brightness temperature and radar backscatter (Larue et al., 2018; Shrestha and Barros, 2025). This critical need is driving improvements in snow representation within LSMs (Navari et al., 2024). In this context, Environment and Climate Change Canada (ECCC), in collaboration with the Canadian Space Agency and Natural Resources Canada, is spearheading the preparation of the Terrestrial Snow Mass Mission (TSMM Derksen et al., 2021). This innovative dual Ku-band spaceborne mission aims to retrieve SWE at a 500 m resolution across the Northern Hemisphere. Currently, ECCC utilizes the Soil Vegetation and Snow version 1.0 (SVS1) (Alavi et al., 2016; Husain et al., 2016; Leonardini et al., 2021) LSM for operational land surface and hydrological predictions across Canada (Gaborit et al., 2017; Durnford et al., 2021). However, SVS1 presents significant limitations: it employs a one-layer snowpack scheme based on the force-restore approach (Leonardini et al., 2021) that does not simulate snow microstructure by design. Moreover, SVS1 does not include key processes affecting snow evolution in forested terrain, such as snow interception by the canopy. These substantial shortcomings underscore the need to enhance SVS1 and significantly expand its capabilities to accurately simulate snow microstructure.

This article presents version 2.0 of SVS (SVS2), focusing specifically on snowpack processes. Built upon the existing code architecture of SVS1, SVS2 introduces major updates across several key components. The snowpack scheme now relies on the externalized version of the detailed Crocus snowpack model (Brun et al., 1992; Vionnet et al., 2012; Lafaysse et al., 2017), incorporating significant developments such as a revised snow albedo scheme (Gaillard et al., 2025), and novel physical parameterizations for Arctic snowpack properties (Woolley et al., 2024) synthesizing earlier works (Barrere et al., 2017; Gouttevin et al., 2018; Royer et al., 2021b). A notable enhancement compared to SVS1 is the new forest canopy module, which includes snow interception and is integrated with Crocus to improve snowpack simulations in forested terrain. Additionally, SVS2 features a new soil module that solves heat diffusion, accounting for soil freezing and thawing processes based on Boone et al. (2000). While SVS2 has been previously applied in several studies (Vionnet et al., 2022; Leroux et al., 2023; Woolley et al., 2024, 2025; Gaillard et al., 2025; Meloche et al., 2024), this article provides a comprehensive description of SVS2 and demonstrates its capability to simulate snowpack properties, including microstructure, across extensive geographical regions. Section 2 gives a general description of SVS2, its snow component and its canopy module. Additional details on non-snow components can be found in the Supplement. Section 3 then details simulation configurations, evaluation datasets, and methodologies. Simulations results are presented and discussed in Sect. 4. Finally, Sect. 5 summarizes the findings and provides concluding remarks.

2 Model description

2.1 Land surface tiling

SVS2 employs a tiling approach to capture land surface heterogeneity within individual grid cells. The model represents six distinct subgrid land surface fractions (“tiles”): (i) exposed snow-free bare ground, (ii) snow-free low vegetation, (iii) snowpack covering bare ground and low vegetation, (iv) high vegetation with intercepted snow, (v) snow-free ground beneath high vegetation canopy, and (vi) snowpack beneath high vegetation canopy. High vegetation in SVS2 encompasses forests of varying density with trees exceeding 3 m in height. Depending on environmental conditions, certain tiles may become inactive within a grid cell (e.g., during snow-free periods or in areas lacking high vegetation). SVS2 computes the surface energy budget and associated surface temperature evolution for each active land surface tile. Figure 1 illustrates the land surface tiling structure in SVS2.

All land surface tiles interact directly with the atmosphere and are energetically and hydrologically coupled with one single atmospheric layer above (Fig. 1). The energy and mass balances of tiles beneath high vegetation (tiles v and vi) experience additional modifications from canopy processes that alter incoming radiation, wind speed, and precipitation through interception mechanisms, as described in Sect. 2.3. Following the approach used in the EC-Land model (Boussetta et al., 2021), SVS2 utilizes a single soil column beneath all land surface tiles. This soil column employs a multi-layer scheme to resolve thermal and hydrological processes, including freeze-thaw dynamics (Boone et al., 2000; Decharme et al., 2011; Alavi et al., 2016). Direct energetic and hydrological coupling exists between the soil column and all surface tiles except high vegetation. High vegetation instead exerts indirect influence through two pathways: (i) modifications to the surface energy budget of underlying tiles (tiles v and vi) and (ii) transpiration that alters soil water content in soil layers.

Snow-covered land surface tiles and the high vegetation tile are described in Sects. 2.2 and 2.3. The rest of SVS2 is described in the Supplement. Section S1.2 in the Supplement details the surface energy budget for exposed and shaded snow-free bare ground and for snow-free low vegetation. Computation of soil thermal regime and hydrology is detailed in Sect. S1.3. Finally, model numerical implementations are described in Sect. S2 and Fig. S2 gives an overview of the order of operations within one model time step.

2.2 Snowpack representation

2.2.1 Detailed snowpack scheme Crocus

In SVS2, snowpack evolution over bare ground, low vegetation, and below high vegetation is simulated with the detailed snowpack scheme Crocus (Brun et al., 1992; Vionnet et al., 2012; Lafaysse et al., 2017). When a SVS2 grid cell contains both high vegetation and open terrain, the Crocus snow model runs twice in parallel for that cell: once in open terrain with unperturbed meteorological forcing and once in forested terrain with energy and mass flux modified by the presence of high vegetation as detailed in Sect. 2.3. The implementation of Crocus within SVS2 relies on the standalone version of Crocus that has been recently developed to facilitate the use of Crocus within existing LSMs. Mazzotti et al. (2024) have used this version of Crocus to develop an hyper-resolution forest snow model.

Crocus is a multilayer snowpack scheme that simulates seasonal evolution of physical properties of the snowpack and its vertical layering. Crocus simulates the main internal snowpack processes, including heat diffusion, compaction, liquid water transport, and snow metamorphism. For each snow layer, Crocus simulates evolution of its thickness, density, liquid water content, temperature, age, and snow microstructure, represented by the snow optical diameter and snow grain sphericity. In addition, Crocus uses a historical variable to identify past presence of liquid water or faceted crystals in the layer. Snow type can be diagnosed from the snow microstructure and historical variables (Brun et al., 1992). Crocus within SVS2 includes the multi-physics ensemble modeling framework ESCROC (Ensemble System Crocus, Lafaysse et al., 2017) that can provide estimates of model uncertainty associated with parameterizations of surface and internal snow process.

Snow optical diameter is inversely proportional to snow specific surface area (SSA) that governs the snowpack optical and microwave properties (e.g. Warren, 1982; Picard et al., 2022). Sphericity is a semi-empirical variable that describes the ratio of angular versus rounded shape in a given snow layer (Brun et al., 1992; Carmagnola et al., 2014). Evolution of optical diameter and of sphericity in Crocus is computed using metamorphism laws as described in Brun et al. (1992), Carmagnola et al. (2014), and Baron (2023). For dry snow in a given layer, the temporal evolution of optical diameter and sphericity are a function of the vertical temperature gradient, whereas, for wet snow, the increase in optical diameter and sphericity with time depends on liquid water content.

Crocus uses a dynamic discretization to represent the evolution of snowpack layering. New snow layers can be created by snowfall events until a user-defined maximum number of layers is reached. Aggregation routines are then used to merge layers with similar physical properties, maintaining sufficiency thin layers close to the snow surface and the bottom of the snowpack to compute heat exchanges with the atmosphere and soil, respectively. A maximum number of 50 snow layers is typically used in the context of operational avalanche hazard forecasting (Vionnet et al., 2012). In the current version of SVS2, the maximum number was fixed at 20 to test a viable configuration for an eventual operational implementation in the context of the Terrestrial Snow Mass Mission (Derksen et al., 2021); such a configuration needs to balance accuracy and computational efficiency. Similarly, Veyssière et al. (2018) used a maximum number of 20 layers when coupling Crocus with a radiative transfer model to simulate microwave emission modeling.

Additional modifications have been made in Crocus within SVS2 to improve simulation of snowpack physical properties across large spatial extents, in particular across the Canadian territory. They include: (i) a revised snow albedo parameterization that takes into account the local climatology of Light-Absorbing Particles (LAP) deposition on snow (Sect. 2.2.2, Gaillard et al., 2025), (ii) additional parameterizations of snow physical processes to simulate properties of Arctic snowpacks (Sect. 2.2.3, Woolley et al. (2024)) and (iii) a canopy module to simulate snow properties in forested environments (Sect. 2.3). Throughout this paper, “SVS2/Crocus” refers specifically to the Crocus version within SVS2, “SURFEX/Crocus” refers to the version within SURFEX, and “Crocus” alone refers to both versions collectively.

2.2.2 Refined snow albedo parametrization

The default snow albedo parameterization in Crocus (Table 1) simulates snow albedo evolution in different large spectral bands: one in the visible (0.3–0.8 µm), and two in the near and shortwave infrared ranges (0.8–1.5 and 1.5–2.8 µm) (Brun et al., 1992; Vionnet et al., 2012). In infrared bands, the spectral albedo depends only on the optical diameter of snow. In the visible band, snow albedo depends on the optical diameter and the amount of LAPs (black carbon, mineral dust) present on the snow surface, which is parameterized by the age of snow and a snow darkening coefficient, γ. This coefficient controls the impact of snow age on the temporal evolution of albedo in the visible band. The default value of γ has been set to 60 d during early developments of Crocus at the Col de Porte experimental site in the French Alps (Brun et al., 1992). Applications of the model to other climates such as Antarctica or the Canadian Arctic, where LAP deposition is significantly lower than in the French Alps, have required calibration of γ (Brun et al., 2011; Woolley et al., 2024).

Table 1 Equations representing the snow albedo for the three spectral bands in the default albedo parameterization of Crocus used in SVS2. The parameters are as follows: d_opt (m) is the optical diameter of the snow, ρ is the snow density, A (days) is the age of the snow, and γ (days) is the snow darkening coefficient. Adapted from Table 4 in Vionnet et al. (2012). The pressure term present when computing Δα_age in SURFEX/Crocus (Vionnet et al., 2012) is not used in SVS2/Crocus as explained in Gaillard et al. (2025).

A local calibration of γ is not possible when applying the model over large domains where deposition of LAP on snow can vary significantly. To overcome this limitation, Gaillard et al. (2025) proposed a new approach for SVS2/Crocus where γ varies in space based on a global climatology of LAP deposition over snow that they developed. Optimized values of γ for multi-year snow albedo simulations with SVS2/Crocus were generated at 10 reference experimental sites spanning a large variety of climates across the world. A regression was then established between these optimal values and the climatological deposition of LAP on snow at the location of the experimental sites. This regression was finally combined with the global climatology to obtain a LAP-informed and spatially variable darkening coefficient for the Crocus albedo parameterization. This methodology of intermediate complexity, implemented in SVS2, considers the spatial variability of LAP deposition in the default snow albedo parametrization without additional computational cost. More details can be found in Gaillard et al. (2025). This updated parameterization improved snow albedo simulations by 10 % at the 10 experimental sites with the largest improvements found in the Arctic (more than 25 %) (Gaillard et al., 2025). The current study represents the first assessment of how this new albedo parameterization affects spatially distributed snowpack simulations (Sect. 4).

Another option for snow albedo calculation, only available in SURFEX/Crocus, relies on the Two-stream Analytical Radiative TransfEr in Snow (TARTES) radiative transfer scheme (Picard and Libois, 2024) and accounts explicitly for the impact of LAP as detailed in Tuzet et al. (2017). This second option is not yet activated in SVS2/Crocus due to a significant additional computational cost (by a factor around 20 for the snow model) and the need for additional atmospheric forcings (LAP deposition fluxes).

2.2.3 Simulation of Arctic snow properties

Crocus was initially developed for avalanche hazard forecasting in alpine terrain (Brun et al., 1992; Durand et al., 1993). Several studies showed that the default version of Crocus does not simulate snowpack properties in Arctic and Subarctic environments well (Domine et al., 2016a; Barrere et al., 2017; Domine et al., 2019; Royer et al., 2021b; Lackner et al., 2022; Woolley et al., 2024), where the Arctic is defined as the area of the Earth with tundra vegetation with the treeline defining the southern limit (Raynolds et al., 2019). In this region, the model fails at simulating the typical vertical profile of snow density (dense wind-slab layer overlaying a layer of lower density made of depth hoar), which negatively impacts simulation of the snow and soil thermal regime in these regions (Barrere et al., 2017; Domine et al., 2019). Using a multi-physics modelling approach, Woolley et al. (2024) have shown that all the physical options available in the default ESCROC ensemble (Lafaysse et al., 2017) are not sufficient to capture the vertical profile of snow density measured at the experimental site of Trail Valley Creek (Northwest Territories, Canada). To overcome this limitation, Woolley et al. (2024) implemented a new set of options in SVS2 to improve the simulation of Arctic snowpack properties based on previous studies (Domine et al., 2016a; Barrere et al., 2017; Gouttevin et al., 2018; Royer et al., 2021b; Lackner et al., 2022). These modifications are referred to as the Arctic version of SVS2/Crocus and described in detail in Woolley et al. (2024). A brief overview is given in the next paragraphs.

The first set of modifications made by Woolley et al. (2024) aims at increasing the compaction of surface snow due to high wind speeds. New parameterizations of snowfall density have been added following Royer et al. (2021b) and Lackner et al. (2022). These parameterizations modify the default parameterization of Vionnet et al. (2012) by increasing the effect of wind speed on falling snow density. Woolley et al. (2024) also adapted the snowdrift scheme of Crocus that simulates wind-packing, the increase in near-surface snow density during drifting and blowing snow events with or without concurrent snowfall (Brun et al., 1997; Vionnet et al., 2012, 2013). This parameterization represents the effect of surface snow fragmentation during wind-induced snow transport and the associated increase in surface snow density (Comola et al., 2017; Sommer et al., 2018). This parameterization has been initially tested in alpine terrain (Vionnet et al., 2013) but evaluations in the Arctic revealed a systematic underestimation of surface density simulated with this approach (Barrere et al., 2017; Royer et al., 2021b; Lackner et al., 2022). Changes to the snowdrift scheme in SVS2/Crocus (Woolley et al., 2024) are based on Royer et al. (2021b) and increase the maximum density of snow impacted by wind and wind-packing rate. When activated, these modifications that increase the compaction of surface snow due to high wind speeds affect simulations over the whole domain (even in regions located outside the Arctic).

The second set of modifications proposed by Woolley et al. (2024) concerns the impact of Arctic vegetation (sedges, shrubs) at the base of the snowpack. In the Arctic version of SVS2/Crocus, snowdrift and blowing snow sublimation schemes are not activated when the total snow height is below the height of Arctic vegetation to simulate that shrubs and other basal vegetation trap wind-blown snow, which limits snow erosion and compaction by wind (Liston et al., 2002; Domine et al., 2016a). The rate of snow compaction in SVS2/Crocus is also reduced below the height of Arctic vegetation (Domine et al., 2016a; Lackner et al., 2022). These two modifications lead to simulation of low density basal snow layers in the presence of Arctic vegetation (Gouttevin et al., 2018; Royer et al., 2021b; Woolley et al., 2024). However, vertical water vapor transport through the snowpack, which contributes to the formation of low-density basal layers in the Arctic (Domine et al., 2019), is not currently included in Crocus. Additional effects of shrubs on winter snow albedo, shortwave vertical irradiance in snow, and the effects of shrub branches on snow metamorphism (Belke-Brea et al., 2020; Domine et al., 2025) are also not simulated by Crocus.

When implementing the effect of Arctic vegetation height in distributed snowpack simulations, Royer et al. (2021b) developed a global map of Arctic vegetation height at 1 km  resolution based on the ECOCLIMAP vegetation database (Masson et al., 2003) used in the SURFEX platform (Masson et al., 2013). The Arctic vegetation database of Royer et al. (2021b) includes alpine-tundra regions located for example in the Rocky Mountains in Canada and the United States. Using this vegetation database, effects of basal vegetation are activated in these alpine-tundra regions. To restrict these effects to the Arctic, a new database of Arctic vegetation height has been developed in SVS2 based on the Circumpolar Arctic Vegetation Map (CAVM Raynolds et al., 2019) at 1 km resolution. This new database is described in Sect. S1.1.3. The CAVM covers the Arctic, above the tree line. Elsewhere, the height of Arctic vegetation is set to zero so that modifications of Woolley et al. (2024) associated with presence of basal vegetation are only active above the Arctic treeline.

2.3 Canopy module

SURFEX/Crocus is coupled to the multiple energy budget scheme (MEB, Boone et al., 2017). MEB-Crocus has been used to simulate snow in forested environments (Nousu et al., 2024). However, MEB is strongly tied to the internal SURFEX code so it is not available in the standalone version of Crocus used in SVS2. For this reason and to address limitations of SVS1 in snow simulations in forests (Leonardini et al., 2021), a new canopy scheme was developed for SVS2. The canopy scheme relies on well-established parameterizations (Hedstrom and Pomeroy, 1998; Koivusalo and Kokkonen, 2002; Gouttevin et al., 2015; Mazzotti et al., 2024; Essery et al., 2025) to represent evolution of the canopy energy and mass balance and its sub-canopy influence on shaded ground and/or snow.

2.3.1 Canopy characteristics

High vegetation in SVS2 refers to forest of varying density covered by tree taller than 3 m. This vegetation is described by three canopy structure parameters: (i) canopy height (h_top), (ii) canopy closure (V_dens) and (iii) average vegetation area index (V_AI) of trees composing the high vegetation. Canopy closure, V_dens, varies between a minimum value of 0.2 (sparse canopy with a canopy coverage of 20 %) and 1 (100 % closed canopy, no gaps). V_AI includes leaves, stems, and branches, and is expressed as one-sided area of woody components and leaves per unit ground of closed forest. The product V_densV_AI converts V_AI per unit ground of closed forest to an effective V_AI per unit ground of the high vegetation tile that accounts for canopy closure. The skyview factor of the canopy (SVF) is derived from V_dens and V_AI (Essery et al., 2008):

$$\begin{array}{}\text{(1)}& \text{SVF}=\text{exp}(-V_{\text{dens}}\mathrm{\Omega }{V}_{\text{AI}})\end{array}$$

where Ω is a clumping parameter that accounts for leaves being vertically stacked in a canopy (Essery et al., 2008; Musselman et al., 2015). In the current version of SVS2, Ω is taken equal to 0.5 as proposed by Bartlett et al. (2006) for needle-leaf trees.

When applied in point-scale mode, h_top and V_dens are provided by the model user and V_AI depends on the vegetation type selected by the user (Table S2). For distributed simulations within GEM-Hydro (Gaborit et al., 2017; Vionnet et al., 2020), h_top is obtained from the 1 km global dataset of Simard et al. (2011) whereas V_dens is derived from the Global Tree Cover dataset (Hansen et al., 2013) at 30 m grid spacing following the methods described in Sect. S1.1.2. Finally, V_AI depends on vegetation types present in the land cover database used as an input for GEM Hydro (Table S2). α_vh and ϵ_vh, the albedo and emissivity of snow-free high vegetation, respectively, also depend on vegetation types as detailed in Table S2.

2.3.2 Canopy energy balance

In SVS2, a single canopy layer that includes leaves, woody parts, and intercepted snow is considered, similar to the one-layer module available in the FSM2 model (Mazzotti et al., 2020; Essery et al., 2025). FSM2 is a reference snowpack model that contains multiple options for representing physical processes in forest canopies. The canopy temperature (T_v, K) in the one-layer module can be determined via:

$$\begin{array}{}\text{(2)}& C_{\mathrm{v}}{\displaystyle \frac{\partial T_{\mathrm{v}}}{\partial t}}=R_{\text{nv}}-H_{\mathrm{v}}-{\text{LE}}_{\mathrm{v}}+L_{\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{v}}\end{array}$$

where C_v is heat capacity of the canopy layer ($\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$) (Eq. 3), R_nv is net radiation budget of the canopy (W m⁻²), H_v and LE_v are sensible and latent heat fluxes from the canopy toward the upper atmosphere (W m⁻²), respectively, L_f is the latent heat of fusion (J kg⁻¹), and Φ_v is the phase change energy of intercepted water on the canopy ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) (Eq. 15). Details about the numerical solving of Eq. (2) are given in Sect. 2.3.

The heat capacity of the canopy, C_v, is written as:

$$\begin{array}{}\text{(3)}& C_{\mathrm{v}}=C_{\text{vb}}+C_{\text{pi}}{S}_{\text{nv}}\end{array}$$

with S_nv the mass of intercepted snow (kg m⁻²), C_pi is specific heat capacity of ice ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) and C_vb = 3.6 × 10⁶ V_AI is the canopy heat capacity ($\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$) (Mazzotti et al., 2020; Essery et al., 2025).

Net radiation to the canopy is the sum of net shortwave radiation (SW_net, W m⁻²) and net longwave radiation (LW_net, W m⁻²). From Gouttevin et al. (2015), net shortwave radiation absorbed by the canopy includes multiple reflections between the canopy and the surface below the canopy:

$$\begin{array}{}\text{(4)}& {\text{SW}}_{\text{net}}={\text{SW}}_{\text{in}}(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{c}}){\mathit{\sigma }}_{\mathrm{f}}\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{s}}(\mathrm{1}-{\mathit{\sigma }}_{\mathrm{f}})}{\mathrm{1}-{\mathit{\alpha }}_{\mathrm{s}}{\mathit{\alpha }}_{\mathrm{c}}{\mathit{\sigma }}_{\mathrm{f}}}}\right)\end{array}$$

with SW_in the incoming solar radiation above the canopy (W m⁻²), α_c the canopy albedo, α_s the albedo of the surface below the canopy, and σ_f the absorption factor (${\mathit{\sigma }}_{\mathrm{f}}=\mathrm{1}-\text{SVF}$). If both bare ground and snow coexist below the canopy, α_s is a weighted average of the snow albedo and the albedo of the bare ground using their fraction areas. The canopy albedo accounts for the presence of intercepted snow:

$$\begin{array}{}\text{(5)}& {\mathit{\alpha }}_{\mathrm{c}}=p_{\text{nv}}{\mathit{\alpha }}_{\text{nv}}+(\mathrm{1}-p_{\text{nv}}){\mathit{\alpha }}_{\text{vh}}\end{array}$$

where α_nv is the albedo of intercepted snow (set to 0.3 as in Gouttevin et al., 2015) and α_vh is the albedo of snow-free high vegetation as defined in Sect. 2.3.1. p_nv is the fraction of canopy covered by intercepted snow (Essery et al., 2003):

$$\begin{array}{}\text{(6)}& p_{\text{nv}}={\left({\displaystyle \frac{S_{\text{nv}}}{S_{\text{nv}}^{\text{max}}}}\right)}^{\mathrm{2}/\mathrm{3}}\end{array}$$

where $S_{\text{nv}}^{\text{max}}$ is the maximum SWE of intercepted snow by canopy over the high vegetation tile (kg m⁻²) as defined in Sect. 2.3.4.

Since emissivities of snow, vegetation and soil are very close to unity, the net longwave radiation budget of the canopy can be simplified as in Gouttevin et al. (2015) and Essery et al. (2025). It is written:

$$\begin{array}{}\text{(7)}& {\text{LW}}_{\text{net}}={\mathit{\sigma }}_{\mathrm{f}}\left({\text{LW}}_{\text{in}}+{\text{LW}}_{\text{surf}}-\mathrm{2}\mathit{\sigma }{T}_{\mathrm{c}}^{\mathrm{4}}\right)\end{array}$$

where LW_in is incoming longwave radiation above the canopy (W m⁻²), σ the Stefan–Boltzman constant, and LW_surf is longwave radiation emitted by the surface below the canopy (W m⁻²). If both bare ground and snow coexist below the canopy, LW_surf is a weighted average of the longwave radiation emitted by snow and longwave radiation emitted by the bare ground using their fraction areas.

The latent and sensible heat fluxes (LE_v and H_v) are computed between the reference level of meteorological inputs (above the canopy) and the canopy surface as in Gouttevin et al. (2015). As a consequence, the latent heat flux between the canopy and upper atmosphere is determined as follows:

$$\begin{array}{}\text{(8)}& {\text{LE}}_{\mathrm{v}}={\mathit{\rho }}_{\mathrm{a}}\left[h_{\mathrm{v}}\right(\mathrm{1}-p_{\text{nv}})L_{\mathrm{v}}+p_{\text{nv}}{L}_{\mathrm{s}}]{\displaystyle \frac{q_{\text{sat\hspace{0.17em}\hspace{0.17em}v}}-q_{\text{air}}}{R_{\text{av}}}}\end{array}$$

where h_v is the Halstead coefficient that partitions total snow-free plant evapotranspiration between evaporation of intercepted water and transpiration via stomata (e.g., Boone et al., 2017) (Eq. 9), ρ_a is the air density (kg m⁻³), R_av is aerodynamic resistance between the canopy and atmosphere above the canopy (s m⁻¹) as derived by the surface turbulent scheme of the GEM model (Beljaars and Holtslag, 1991; Delage and Girard, 1992), q_sat  v is specific humidity at saturation of the canopy layer (kg kg⁻¹), q_air is specific humidity of the atmosphere above the canopy (kg kg⁻¹) and L_v and L_s are the latent heat of vaporization and sublimation (J kg⁻¹), respectively.

The Halstead coefficient, h_v, is computed as:

$$\begin{array}{}\text{(9)}& h_{\mathrm{v}}=\left({\displaystyle \frac{R_{\text{av}}}{R_{\text{av}}+R_{\mathrm{s}}}}\right)(\mathrm{1}-\mathit{\delta })+\mathit{\delta }\end{array}$$

where R_s is the stomatal resistance (s m⁻¹) obtained from the SVS photosynthesis module (Husain et al., 2016) and δ is the fraction of intercepted liquid water (Deardorff, 1978) computed as:

$$\begin{array}{}\text{(10)}& \mathit{\delta }={\left({\displaystyle \frac{W_{\text{wv}}}{W_{\text{wv}}^{\text{max}}}}\right)}^{\mathrm{2}/\mathrm{3}}\end{array}$$

where W_wv and $W_{\text{wv}}^{\text{max}}$ are intercepted liquid water and maximum capacity of intercepted liquid water in the canopy over the high vegetation tile (kg m⁻²), respectively (Sect. 2.3.4).

Total water flux from the canopy to the atmosphere (E_v, $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, Eq. 11) can be decomposed into the transpiration rate of vegetation (E_tr, $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, Eq. 12), the evaporation rate of intercepted liquid water (E_r, $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, Eq. 13), and sublimation of intercepted snow (E_rs, $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, Eq. 14).

$$\begin{array}{}\text{(11)}& E_{\mathrm{v}}={\text{LE}}_{\mathrm{v}}/L=E_{\text{tr}}+E_{\mathrm{r}}+E_{\text{rs}}\end{array}$$

where L is the effective latent heat (J kg⁻¹) as in Boone et al. (2017) ($L=p_{\text{nv}}{L}_{\mathrm{s}}+(\mathrm{1}-p_{\text{nv}})L_{\mathrm{v}}$). E_tr, E_r and E_rs are written as:

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }{E}_{\text{tr}}={\mathit{\rho }}_{\mathrm{a}}(L_{\mathrm{v}}/L){\displaystyle \frac{q_{\text{sat\hspace{0.17em}\hspace{0.17em}v}}-q_{\text{air}}}{R_{\text{av}}+R_{\mathrm{s}}}}(\mathrm{1}-\mathit{\delta })(\mathrm{1}-p_{\text{nv}})\text{(13)}& {\displaystyle }& {\displaystyle }{E}_{\mathrm{r}}={\mathit{\rho }}_{\mathrm{a}}(L_{\mathrm{v}}/L){\displaystyle \frac{q_{\text{sat\hspace{0.17em}\hspace{0.17em}v}}-q_{\text{air}}}{R_{\text{av}}}}\mathit{\delta }(\mathrm{1}-p_{\text{nv}})\text{(14)}& {\displaystyle }& {\displaystyle }{E}_{\text{rs}}={\mathit{\rho }}_{\mathrm{a}}(L_{\mathrm{s}}/L){\displaystyle \frac{q_{\text{sat\hspace{0.17em}\hspace{0.17em}v}}-q_{\text{air}}}{R_{\text{av}}}}{p}_{\text{nv}}.\end{array}$$

The term of phase change energy due to melt of intercepted snow or refreezing of intercepted water in canopy is determined following Boone et al. (2017):

$$\begin{array}{}\text{(15)}& {\mathrm{\Phi }}_{\mathrm{v}}=S_{\text{nv}}{k}_{{\mathrm{\Phi }}_{\mathrm{v}}}(T_{\mathrm{m}}-T_{\mathrm{v}})\end{array}$$

where $k_{{\mathrm{\Phi }}_{\mathrm{v}}}$ is a constant equal to 5.56 × 10⁻⁶ ${\mathrm{s}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (Boone et al., 2017) and T_m is the freezing temperature of liquid water (K).

Sensible heat flux between the canopy and atmosphere is determined as followed:

$$\begin{array}{}\text{(16)}& H_{\mathrm{v}}={\mathit{\rho }}_{\mathrm{a}}{c}_p{\displaystyle \frac{T_{\mathrm{v}}-T_{\mathrm{a}}}{R_{\text{av}}}}\end{array}$$

where c_p is specific heat of dry air ($\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{kg}}^{-\mathrm{1}}$).

2.3.3 Wind profile and aerodynamic resistances

The canopy modifies the wind profile close to the surface and affects turbulent exchanges between the canopy and atmosphere, and the surface below the canopy and atmosphere. The wind profile is assumed to be logarithmic above the canopy (including stability correction), exponential within the canopy (no stability correction), and again logarithmic below the canopy (no stability correction) as often done in land surface and snow models (e.g., Koivusalo and Kokkonen, 2002; Mahat et al., 2013; Mazzotti et al., 2024; Essery et al., 2025):

$$\begin{array}{}\text{(17)}& U\left(z\right)=\left\{\begin{array}{ll}{U}_{\mathrm{a}}\left[\text{ln}\left(\frac{z-d}{z_{\mathrm{0}\text{mv}}}\right)-{\mathrm{\Psi }}_{\mathrm{m}}\left(\frac{z-d}{L_{\text{MO}}}\right)\right.& \\ \left.+{\mathrm{\Psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{0}\text{mv}}}{L_{\text{MO}}}\right)\right]\left[\text{ln}\left(\frac{z_{\mathrm{U}}-d}{z_{\mathrm{0}\text{mv}}}\right)\right.& \\ -{\mathrm{\Psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{U}}-d}{L_{\text{MO}}}\right)& \\ {\left.+{\mathrm{\Psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{0}\text{mv}}}{L_{\text{MO}}}\right)\right]}^{-\mathrm{1}},& \text{if }z\ge h_{\text{top}}\\ U_{\text{top}}\text{exp}\left[w_{\text{can}}\left(\frac{z}{h_{\text{top}}}-\mathrm{1}\right)\right],& \text{if }{h}_{\text{bot}}<z\\ & <h_{\text{top}}\\ U_{\text{bot}}\text{ln}\left(\frac{z}{z_{\mathrm{0}}}\right){\left[\text{ln}\left(\frac{h_{\text{bot}}}{z_{\mathrm{0}}}\right)\right]}^{-\mathrm{1}},& \text{if }z\le h_{\text{bot}}\end{array}\right.\end{array}$$

where h_top is canopy height (m), h_bot is canopy base height (set constant at 2 m), U_a is wind speed (m s⁻¹) at the forcing height z_U (m) (above the ground, z_U>h_top), U_top and U_bot are wind speed at the top and bottom of the canopy (m s⁻¹), z_0mv is momentum roughness length of the canopy (m) and d is the zero plane displacement height (m). z_0mv and d depend on the canopy height following: z_0mv = 0.1 h_top and d = 0.67 h_top. The factor w_can characterizes attenuation of wind speed in the canopy and depends on canopy structure with w_can=βV_densV_AI, with β = 0.9 as in Liston and Elder (2006). The Monin–Obukhov length, L_MO, and the stability function for momentum, Ψ_m, are computed by the surface turbulence scheme of the GEM atmospheric model (Beljaars and Holtslag, 1991; Delage and Girard, 1992).

Aerodynamic resistance between the canopy and atmosphere (reference level above the canopy), R_av and resistance between the surface below the canopy and the atmosphere (reference level above the canopy), R_bv, are computed using the diffusivity theory (Koivusalo and Kokkonen, 2002; Mahat et al., 2013; Boone et al., 2017; Essery et al., 2025). In the canopy, eddy diffusivity, K_c, follows the same exponential form as the wind profile:

$$\begin{array}{}\text{(18)}& K_{\mathrm{c}}=K\left(h_{\text{top}}\right)\text{exp}\left[w_{\text{can}}\left({\displaystyle \frac{z}{h_{\text{top}}}}-\mathrm{1}\right)\right]\end{array}$$

where K(h_top) is the eddy diffusivity at top of the canopy (m² s⁻¹) and is defined as:

$$\begin{array}{}\text{(19)}& K\left(h_{\text{top}}\right)=\mathit{\kappa }{u}_{{}^{\ast }\text{top}}(h_{\text{top}}-d)\end{array}$$

where κ is the von Karman constant (–) and $u_{{}^{\ast }\text{top}}$ is friction velocity at the top of the canopy (m s⁻¹) (including stability correction):

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{u}_{{}^{\ast }\text{top}}=& \mathit{\kappa }{U}_{\mathrm{a}}\left[\text{ln}\left({\displaystyle \frac{z_{\mathrm{U}}-d}{z_{\mathrm{0}\text{mv}}}}\right)-{\mathrm{\Psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{U}}-d}{L_{\text{MO}}}}\right)\right.\\ & {\left.+{\mathrm{\Psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{0}\text{mv}}}{L_{\text{MO}}}}\right)\right]}^{-\mathrm{1}}.\end{array}\end{array}$$

R_av is obtained by integrating the inverse of the eddy diffusion coefficient over the range from d+z_0mv to z_T, the height of atmospheric forcing for air temperature and humidity (Koivusalo and Kokkonen, 2002; Essery et al., 2025):

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{R}_{\text{av}}=& {\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }{u}_{{}^{\ast }\text{top}}}}\left[\text{ln}\left({\displaystyle \frac{z_{\mathrm{T}}-d}{h_{\text{top}}-d}}\right)-{\mathrm{\Psi }}_{\mathrm{H}}\left({\displaystyle \frac{z_{\mathrm{T}}-d}{L_{\text{MO}}}}\right)\right.\\ & \left.+{\mathrm{\Psi }}_{\mathrm{H}}\left({\displaystyle \frac{h_{\text{top}}-d}{L_{\text{MO}}}}\right)\right]+{\displaystyle \frac{h_{\text{top}}}{w_{\text{can}}K\left(h_{\text{top}}\right)}}\\ & \left[\text{exp}\left(w_{\text{can}}-{\displaystyle \frac{w_{\text{can}}(d+z_{\mathrm{0}\text{mv}})}{h_{\text{top}}}}\right)-\mathrm{1}\right]\end{array}\end{array}$$

Aerodynamic resistance between the surface below the canopy and the atmosphere, R_bv, is calculated as in Koivusalo and Kokkonen (2002) and Essery et al. (2025):

$$\begin{array}{}\text{(22)}& R_{\text{bv}}=R_{\text{av}}+R_{\text{surf}}\end{array}$$

where R_surf combines resistance of the lower part of the canopy (h_bot to d+z_0mv) and resistance between the surface and the canopy bottom (z₀ to h_bot):

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{R}_{\text{surf}}=& {\displaystyle \frac{h_{\text{top}}{e}^{w_{\text{can}}}}{w_{\text{can}}K\left(h_{\text{top}}\right)}}\left[\text{exp}\left(-{\displaystyle \frac{w_{\text{can}}{h}_{\text{bot}}}{h_{\text{top}}}}\right)\right.\\ & \left.-\text{exp}\left(-{\displaystyle \frac{w_{\text{can}}(d+z_{\mathrm{0}\text{mv}})}{h_{\text{top}}}}\right)\right]\\ & +{\displaystyle \frac{\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}{U}_{\text{bot}}}}\text{ln}\left({\displaystyle \frac{h_{\text{bot}}}{z_{\mathrm{0}\mathrm{h}}}}\right)\text{ln}\left({\displaystyle \frac{h_{\text{bot}}}{z_{\mathrm{0}}}}\right)\end{array}\end{array}$$

where z₀ and z_0h are roughness lengths for momentum and heat for the surface below the vegetation (shaded ground or snow).

2.3.4 Mass balance of snow and liquid water in the canopy

Liquid water (W_wv, kg m⁻²) and snow (S_nv, kg m⁻²) can exist on the canopy. Both quantities are calculated with respect to the total high vegetation tile area. The mass balance of intercepted liquid water follows:

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\partial W_{\text{wv}}}{\partial t}}=(\mathrm{1}-\text{SVF})\text{PR}-V_{\text{dens}}{E}_{\mathrm{r}}-D_{\mathrm{r}}-{\mathrm{\Phi }}_{\mathrm{v}}\end{array}$$

where PR is the rate of liquid precipitation falling above the canopy ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) and SVF is the sky view factor (Eq. 1), quantifying how much of the precipitation falls directly onto the ground. Φ_v is the source/sink term associated with melt of intercepted snow or refreezing of intercepted water (Eq. 15, $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). A positive value of Φ_v corresponds to freezing of intercepted water. Finally, D_r is the amount of liquid water exceeding the maximum capacity of intercepted liquid water in the canopy, $W_{\text{wv}}^{\text{max}}$, that drips below the canopy ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). Numerically, D_r is computed as:

$$\begin{array}{}\text{(25)}& D_{\mathrm{r}}=\text{max}\left(\mathrm{0},W_{\text{wv}}-W_{\text{wv}}^{\text{max}}\right)/\mathrm{\Delta }t\end{array}$$

where $W_{\text{wv}}^{\text{max}}=\mathrm{0.2}{V}_{\text{dens}}{V}_{\text{AI}}$ as in Dickinson (1984) and Husain et al. (2016),

The rate of liquid precipitation falling on the surface below the canopy ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) includes dripping liquid water and rain falling through the canopy:

$$\begin{array}{}\text{(26)}& {\text{PR}}_{\mathrm{c}}=\text{SVF}\times \text{PR}+D_{\mathrm{r}}.\end{array}$$

The mass balance of intercepted snow in the canopy follows:

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\partial S_{\text{nv}}}{\partial t}}=I_{\text{sv}}-U_{\text{sv}}-V_{\text{dens}}{E}_{\text{rs}}+{\mathrm{\Phi }}_{\mathrm{v}}\end{array}$$

where I_sv is the interception rate of snowfall by the canopy (Eq. 28, $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), E_rs the sublimation rate of intercepted snow (Eq. 14, $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), Φ_v is the source/sink term associated with the melt of intercepted snow or refreezing of intercepted water (Eq. 15, $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) and U_sv is the snow unloading rate from the canopy ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). U_sv follows the simple approach of Mazzotti et al. (2020) with $U_{\text{sv}}=S_{\text{nv}}\frac{\mathrm{\Delta }t}{\mathit{\tau }}$, with Δt the model time step (s) and τ a relaxation time constant (s). This time constant varies with the canopy temperature such that when the canopy temperature is below freezing, this time constant is taken equal to 240 h, otherwise it equals to 48 h.

The amount of snowfall intercepted by the canopy, I_sv, is taken from Hedstrom and Pomeroy (1998):

$$\begin{array}{}\text{(28)}& I_{\text{sv}}=\left(S_{\text{nv}}^{\text{max}}-S_{\text{nv}}\right)/\mathrm{\Delta }t\left(\mathrm{1}-\text{exp}\left(-V_{\text{dens}}\text{SR}\mathrm{\Delta }t/W_{\text{sv}}^{\text{max}}\right)\right)\end{array}$$

where SR is the snowfall rate ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) and $S_{\text{nv}}^{\text{max}}$ is the maximum snow holding capacity of the canopy. $S_{\text{nv}}^{\text{max}}$ varies with the air temperature (T_a, °C) as in Andreadis et al. (2009):

$$\begin{array}{}\text{(29)}& S_{\text{nv}}^{\text{max}}=\left\{\begin{array}{ll}\mathrm{4}{m}_{\text{scap}}{V}_{\text{dens}}\text{LAI},& \text{if }{T}_{\mathrm{a}}>-\mathrm{1}\mathrm{°}\mathrm{C}\\ (\mathrm{1.5}{T}_{\mathrm{a}}-\mathrm{5.5})m_{\text{scap}}& \\ \times V_{\text{dens}}\text{LAI},& \text{if }-\mathrm{3}\mathrm{°}\mathrm{C}<T_{\mathrm{a}}\\ & \le -\mathrm{1}\mathrm{°}\mathrm{C}\\ m_{\text{scap}}{V}_{\text{dens}}\text{LAI}& \text{otherwise}\end{array}\right.\end{array}$$

where m_scap is a coefficient fixed to 5 kg m⁻² (Andreadis et al., 2009). The dependency on air temperature of $S_{\text{nv}}^{\text{max}}$ allows the model to capture the increase in snow loading resulting from increased adhesion and cohesion above −3 °C as recommended by Lundquist et al. (2021).

The rate of solid precipitation falling on the surface below the canopy ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) includes the unloading of snow and snowfall not intercepted by the canopy:

$$\begin{array}{}\text{(30)}& {\text{SR}}_{\mathrm{c}}=U_{\text{sv}}+\text{SR}-I_{\text{sv}}\end{array}$$

2.3.5 Impact on sub-canopy energy and mass fluxes

The presence of canopy in SVS2 modifies the energy and mass fluxes received by the land surface tiles below the canopy (shaded bare ground and/or snow). For radiative fluxes, longwave and shortwave radiation received at the surface below the canopy can be expressed as:

$$\begin{array}{}\text{(31)}& {\displaystyle }& {\displaystyle }{\text{LW}}_{\text{bv}}={\text{LW}}_{\text{in}}{\mathit{\tau }}_{\text{dif}}+(\mathrm{1}-{\mathit{\tau }}_{\text{dif}})\mathit{\sigma }{T}_{\mathrm{v}}^{\mathrm{4}}\text{(32)}& {\displaystyle }& {\displaystyle }{\text{SW}}_{\text{bv}}={\text{SW}}_{\text{dir}}{\mathit{\tau }}_{\text{dir}}+{\text{SW}}_{\text{dif}}{\mathit{\tau }}_{\text{dif}}\end{array}$$

where LW_in is longwave radiation above the canopy (W m⁻²), SW_dir and SW_dif are direct and diffuse shortwave radiation above the canopy (W m⁻²), T_v the canopy temperature (K) of the previous time step due to the sequential approach used in SVS2, τ_dif is the transmissivity of diffuse radiation through the canopy, which is taken equal to the skyview factor, SVF (Eq. 1), as in Musselman et al. (2015), and τ_dir is the direct transmissivity of shortwave radiation defined as in Pomeroy and Dion (1996):

$$\begin{array}{}\text{(33)}& {\mathit{\tau }}_{\text{dir}}={\text{e}}^{-k_{\text{ext}}{V}_{\text{dens}}{V}_{\text{AI}}\mathrm{\Omega }/\text{sin}\left(\mathit{\theta }\right)}\end{array}$$

where k_ext is an extinction coefficient taken equal to 0.5 for randomly oriented canopy elements and θ is the solar angle (rad). If separation between direct and diffuse shortwave is not available in the meteorological forcing, it is recomputed in the model following Erbs et al. (1982) and Jonas et al. (2020).

Turbulent fluxes of sensible and latent heat above the surface below high vegetation are computed using air temperature and specific humidity taken at the reference forcing level above the canopy as in Koivusalo and Kokkonen (2002) and Gouttevin et al. (2015). The canopy effect on turbulent heat fluxes is accounted for in the computation of R_bv, the aerodynamic resistance between the surface below the canopy and the reference forcing level above the canopy (Eq. 22). For the snowpack below high vegetation, the wind speed below the canopy, U_bot (Eq. 17), is used in SVS2/Crocus to compute properties of falling snow, wind-induced compaction and mass loss due to blowing snow sublimation. The canopy temperature, T_v, is not used when computing the turbulent fluxes below high vegetation and only affects the incoming longwave radiation below high vegetation (Eq. 31).

Precipitation falling on the surface below the canopy is computed by the canopy mass balance module (Eqs. 26 and 30 in Sect. 2.3.4). Snowfall through the canopy and the snow unloading rate are provided as two distinct input mass fluxes to Crocus. The properties (density, microstructure) of new snow from snowfall are derived from the Crocus parameterization (Vionnet et al., 2012; Lafaysse et al., 2017) whereas the properties of unloaded snow differs from of the properties of snowfall as recommended by Bouchard et al. (2024). SVS2 assumes a fixed density of 200 kg m⁻³, an optical diameter of 0.2 mm and a sphericity of 0.95 for unloaded snow.

2.3.6 Alternatives to the canopy modules

The canopy energy balance and mass balance modules described in Sects. 2.3.2 and 2.3.4 are activated when meteorological forcings are provided at a given height above the canopy. This option is referred to as the ABV option (“above”) in SVS2 and is used when the model is driven by atmospheric forcing taken from meteorological models, reanalysis, or towers located above the canopy. Two other options have been added in SVS2 to allow the model to be driven by atmospheric forcing collected at meteorological stations that are not always located above the canopy. The first option is the “forest” option (referred to as FOR) when the model is driven by meteorological forcing measured below the canopy. In this case, mass and energy processes associated with the canopy are not simulated as all the atmospheric forcings are valid below the canopy and already account for the effect of the canopy. The second option is the “open to forest” option (referred to as O2F) when the model is driven by meteorological forcing transferred from adjacent open terrain. An example would be when the meteorological forcing is collected at 2 m above the ground in an adjacent clearing. This method is often used in snowpack models (Marsh et al., 2020; Mazzotti et al., 2024; Strasser et al., 2024). In SVS2, this approach is only used for model evaluation at point-scale with meteorological forcing data collected in open terrain (e.g. Rutter et al., 2009; Bonner et al., 2022b).

In the O2F option, the canopy energy balance module is not activated and the model assumes that air temperature and relative humidity close to the surface in the open are equal to those below the canopy. The model also assumes that the canopy temperature, T_v is equal to the air temperature in the open (Rutter et al., 2023). Under these assumptions, shortwave and longwave radiation below the canopy are derived using Eqs. (32) and (31), respectively. Wind speed close to the surface in open terrain is assumed equal to wind speed above the canopy and is transferred to the forcing height below the canopy using the wind profile described in Eq. (17). In the O2F option, stability effects above the canopy cannot be computed and the model assumes a logarithmic wind profile above the canopy. Finally, modified meteorological forcing below the canopy are used to drive the energy and mass balance of shaded bare ground and snowpack below high vegetation.

In the O2F option, the canopy mass balance module is activated and simulates snow interception and associated unloading and sublimation as described in Sect. 2.3.4. However, since the energy balance of the canopy is not calculated, the model cannot derive sublimation of intercepted snow from Eq. (14). As a consequence, sublimation of intercepted snow is estimated following Pomeroy et al. (1998) and Essery et al. (2003) as described in Sect. S1.4 of the Supplement. In the O2F option, the source/sink term associated with melt of intercepted snow or refreezing of intercepted water (Φ_v in Eqs. 24 and 27) is not considered and the model relies on simple assumptions similar to Mazzotti et al. (2024). When air temperature drops below freezing, the model treats any liquid water on the canopy as frozen, incorporating it into intercepted snow. Conversely, when air temperature rises above freezing, unloaded snow is treated as liquid water. This meltwater first joins the existing reservoir of intercepted liquid water on the canopy, and any excess beyond the canopy's maximum water-holding capacity will eventually drip to the ground.

3 Methods for model evaluation and benchmarking

3.1 Overview

This section describes two simulation strategies used to evaluate SVS2/Crocus performance. First, point-scale simulations were performed to investigate the canopy module and its effects on simulated snowpack properties. These simulations also included comparisons with SVS1 to demonstrate model improvements. Second, spatially distributed simulations covering a domain spanning from southern Quebec to the Canadian Arctic were conducted to evaluate both the revised albedo parameterization in SVS2/Crocus and the Arctic snow parameterizations. Prior to this work, these modifications had been tested only at the point-scale. The distributed simulations provided an opportunity to assess these parameterizations across large spatial scales and diverse snowpack types within the domain. These simulations enabled us to (i) quantify how the albedo parameterization of Gaillard et al. (2025) affects snow melt dynamics throughout the domain and (ii) evaluate the spatial applicability of both the default and Arctic configurations of SVS2/Crocus for estimating vertical profiles of snow properties (SSA and density). Results from the distributed simulations were also compared against reference snow products used as benchmarks.

3.2 Point-scale simulations

3.2.1 Snow Crested Butte

An evaluation of SVS2 in point-scale mode was conducted at Snow Crested Butte (Colorado, USA) between October 2018 and June 2021. Data used for forcing and evaluating the model were quality controlled and are fully presented in Bonner et al. (2022b). Snow Crested Butte has two adjacent meteorological stations, one in a forest and one in a clearing, about 300 m apart. Data from the clearing (collected at ∼ 2 m above the ground) were used to drive the model in both the open and in the forest using the O2F approach. Following Bonner et al. (2022b), a dew point temperature threshold of 0 °C was used to separate total precipitation into liquid or solid precipitation. Snow pits (vertical profiles of snow density, temperature, and grain type) and snow transects (bulk SWE and snow height) were conducted at both locations and were used to evaluate the simulations. Canopy closure was taken equal to 60 % based on the Global Tree Cover dataset (Hansen et al., 2013). For this test case, SVS2/Crocus relied on the default configuration of the model (Lafaysse et al., 2017). The updated snow albedo parameterization of Gaillard et al. (2025) (Sect. 2.2.2) and changes of Woolley et al. (2024) specific to arctic environments (Sect. 2.2.3) were not activated since no impact was expected in this environment. Simulations were also carried out with SVS1 using the configuration detailed in Leonardini et al. (2021).

3.2.2 Additional tests

Additional tests in point-scale mode focused on specific physical processes affecting the snow cover. SVS2 was tested at the Marmot Creek Research Basin, Kananaskis Country (Alberta, Canada) (Fang et al., 2019) to evaluate the impact of the configuration of the canopy module on (i) the quality of the simulations of bulk snowpack properties (SWE and snow height) in a forest stand compared to simulations in an adjacent clearing, and (ii) the estimation of the radiative fluxes below high vegetation. This evaluation is presented in Sect. S3.1. A direct evaluation of snow interception and unloading is also presented in Sect. S3.2.

Figure 2 Maps of the study domain showing (a) the forest fraction derived from the North American Land Change Monitoring System (NALCMS) dataset and (b) the height of Arctic vegetation. The limits of the simulation domain are shown as thick orange lines on both maps. The locations of the sites where snow pit measurements were collected in 2007/2008 and 2015/2016 are also indicated on both figures.

3.3 Distributed simulations across a latitudinal transect in Eastern Canada

3.3.1 Study domain and period

The ability of SVS2 to simulate large-scale snowpack properties and dynamics was evaluated over a domain that extends from southern Quebec to the North of the Canadian Arctic (43 to 84° N, Fig. 2). This domain corresponds roughly to the latitudinal transect studied by Royer et al. (2021a) when proposing a new classification for northern snowpacks. It includes transitions between several ecozones from boreal forest to subarctic and arctic ecosystems (Fig. 2), and covers three main types of snowpack (Royer et al., 2021a): boreal forest snow (between 47 and 58° N), tundra snow (between 58 and 74° N) and polar desert snow (above 74° N). Two winters have been selected to simulate snowpack evolution over the study domain. The first winter (2007/2008) corresponds to the International Polar Year when snow properties across a large latitudinal spatial extent were measured in Quebec (Langlois et al., 2010; Royer et al., 2021a) (Fig. 2). As no SSA measurements were available in 2007/2008, a second year (2015/2016) is included with vertical profiles of SSA and density for 7 locations across the transect (Fig. 2b) extending from Forêt Montmorency (47.31° N) in southern Quebec up to Ward Hunt (83.09° N) at the northernmost point of the Canadian Arctic (Royer et al., 2021a).

Vionnet et al. (2012)Vionnet et al. (2013)Brun et al. (1992)Vionnet et al. (2012)Vionnet et al. (2013)Brun et al. (1992)Royer et al. (2021b)Royer et al. (2021b)Royer et al. (2021b)

Table 2 Configurations of SVS2 used for the distributed simulations. The LAP-informed values of the darkening coefficient were obtained from Gaillard et al. (2025).

3.3.2 Model forcing and configurations

SVS2 was used in offline (uncoupled) mode to simulate snow cover evolution. For each winter, snowpack simulations were carried out from 1 September to 31 August over the study domain shown on Fig. 2 at 10 km horizontal grid spacing. The grid is made of 210 by 450 points. Soil properties (sand/clay percentages) were obtained from the Global Soil Dataset for use in Earth System Models (GSDE, Shangguan et al., 2014) at 30 arcsec, whereas information about the land cover (vegetation type/fraction) were derived from the NALCMS dataset at 30 m grid spacing (North American Land Change Monitoring System; http://www.cec.org/north-american-land-change-monitoring-system/, last access: 13 November 2025). The influence of soil organic matter on soil thermal and hydraulic properties (Sect. S1.3.4) was not considered. The initial soil state on 1 September of each winter was taken from a 1 year spin-up of the land surface conditions. Such spin-up duration is short compared to typical spin-up durations (e.g., hundreds of years) used in studies focusing on the thermal regime of the permafrost at high-latitudes (e.g., Elshamy et al., 2020; Ji et al., 2022). In the context of this study, a sensitivity analysis showed that results for seasonal snow (snow physical properties and snow melt out date) remained unchanged when using a 5 year spin-up of the land surface conditions.

SVS2 requires the following atmospheric forcing: air temperature and specific humidity at a known level above the surface, wind speed at a known level above the surface, surface air pressure, surface incoming longwave and total shortwave (direct and diffuse) radiations and solid and liquid precipitation rates. These hourly fields were obtained from the ERA5 global reanalysis (Hersbach et al., 2020). Air temperature/specific humidity and wind speed were taken at the lowest prognostic levels (∼ 30 m). Downscaling of the atmospheric forcing (Bernier et al., 2011) was applied to account for different surface elevations between ERA-5 forcing at a horizontal resolution of 31 km and the SVS2 grid at a horizontal resolution of 10 km. The method of Bernier et al. (2011) modifies air temperature of the ERA-5 forcing using a constant lapse rate, adjusts air pressure and assumes that relative humidity remains unchanged. The specific humidity is finally recomputed to account for elevation differences. Precipitation phase was obtained following the method of Harder and Pomeroy (2013), relying on hydrometeor temperature derived from the downscaled air temperature and humidity. The method of Harder and Pomeroy (2013) has been extensively used when simulating cold-region hydrology across Canada (e.g. Aygün et al., 2022; Leroux et al., 2023; Marsh et al., 2024). A similar downscaling approach is used to produce the atmospheric forcing for the ERA5-Land surface reanalysis at 9 km grid spacing (Muñoz-Sabater et al., 2021). The ABV approach was used for high vegetation meaning that the canopy module in SVS2 was fully activated.

Three configurations of SVS2 were considered and are described in Table 2. Two configurations (Def, Def_LAP) were used to first assess the impact of the LAP-informed snow darkening coefficient on snowpack dynamics. An additional configuration (Arc_LAP) was then considered to test the capacity of Arctic configuration (Sect. 2.2.3) to simulate snowpack properties (density and SSA) across the latitudinal transect. Arc_LAP accounts for wind-induced snow compaction through a modified snowfall density and snowdrift routine and includes effects of basal vegetation that depend on the height of Arctic vegetation (Sect. 2.2.3 and Fig. 2). Arc_LAP also includes the LAP-informed snow darkening coefficient. The configuration Def was used for the 1 year spin-up simulations. All simulations include the parameterization of mass loss due to blowing snow sublimation of Gordon et al. (2006) that is recommended for large scale applications of Crocus (Brun et al., 2013). The three configurations use the formulation of Yen (1981) for the snow thermal conductivity.

3.3.3 Evaluation, benchmarking data, and performance metrics

Simulations of internal snowpack properties were evaluated using snow pit data aggregated in Royer et al. (2021a) and collected by different research groups (Langlois et al., 2010; Domine et al., 2018; Lackner et al., 2022). These data consist of vertical profiles of snow density, SSA (for winter 2015/16), snow temperature, and grain type. Only snow density and SSA were considered in this evaluation. A vegetation type (forest or open) was associated with each snow pit. For model evaluation, the simulated profiles of snow density and optical diameter were then extracted at the nearest grid points to each snow pit. Simulated optical diameter, d_opt (m), was converted to SSA (m² kg⁻¹) for the evaluation following:

$$\begin{array}{}\text{(34)}& \text{SSA}={\displaystyle \frac{\mathrm{6}}{d_{\text{opt}}{\mathit{\rho }}_{\text{ice}}}}\end{array}$$

where ρ_ice is the ice density (kg m⁻³). The SVS2 simulated snow pit in the open (forest) was considered when the observation was collected in open (forested) terrain. Profiles where the simulated and observed snow heights differed by more than 50 % were not considered to limit the errors associated with site representativeness. For example, snow profiles collected in deep snow drift in the Arctic were not considered since the model configuration cannot reproduce such specific environments. 32 snow pits out of 203 were excluded using this criterion. For the selected profiles, simulated layer thicknesses were uniformly scaled to match the observed snow depth. Both the observed and scaled simulated profiles were then re-sampled onto a common 1 mm vertical grid to enable direct comparison, following the approach of Woolley et al. (2024). This normalization represents the initial step in the evaluation framework proposed by Lehning et al. (2001) and Viallon-Galinier et al. (2020) for comparing simulated and observed snow profiles. Error metrics (bias and root-mean-square error, RMSE) were computed using normalized profiles for the whole snowpack, the top 25 %, and the bottom 25 %.

The snow melt-out date (SMOD) was also used to assess model performance and compare it against reference benchmarks. For winters 2007/08 and 2015/16, the observed SMOD was estimated using the Global Multisensor Automated Snow and Ice Mapping System (GMASI) product (Romanov, 2017). Based on a combination of satellite observations in the visible, infrared, and microwave bands, GMASI provides daily global estimates of snow presence at 0.04^∘ resolution. SMOD was computed as in Gascoin et al. (2019) and represents the last date of the longest continuous snow period. For the SVS2 experiments (Table 2), a minimum threshold of 5 cm was applied to the grid-averaged snow heights to determine the presence of snow in a grid cell as in Garnaud et al. (2021). Two snow products were used as additional reference benchmarks to evaluate SVS2 simulation performance (i) the ERA5-Land surface reanalysis (Muñoz-Sabater et al., 2021) at 9 km resolution (referred as ERA5L) and (ii) the ERA5-Crocus dataset (Decharme and Barbu, 2024; Ramos Buarque et al., 2025) that covers the Northern Hemisphere at 25 km resolution (referred as ERA5_CRO). Both datasets use the same atmospheric forcing as SVS2 simulations. ERA5L relies on the EC-Land land surface scheme (Boussetta et al., 2021) while ERA5_CRO uses the SURFEX/Crocus (Vionnet et al., 2012) in the configuration described by Brun et al. (2013). ERA5_CRO relies on the same Crocus configuration as Def of SVS2/Crocus (Table 2) and uses 60 d for the snow darkening coefficient (Sect. 2.2.2). The same snow height threshold of 5 cm was applied to the simulations from ERA5L and ERA5_CRO to determine snow presence and derive associated SMOD. Finally, all SMOD data were interpolated to the 10 km grid of SVS2/Crocus.

Figure 3 Simulated (with SVS1 in blue and SVS2 in green) and observed snow height, SWE, and density between October 2019 and June 2021 at (a, c, and e) the open site and (b, d, and f) the forest site at Snow Crested Butte (USA). The black line in (a, b) represents the snow height measured by automatic sensors at each station, the circles represent the mean measurements and the + signs represent the minimum and maximum measurements.

Table 3 Error metrics (Bias and RMSE) for the simulations of bulk snowpack properties (snow depth, SWE and bulk density) with SVS1 and SVS2 at Snow Crested Butte between October 2019 and June 2021. Measurements from automatic stations, snow courses and snow pits have been used to compute the errors metrics. Values in bold represent the best error metric among the two models for each variable and for each location (open and forest).

4 Results and discussions

4.1 Point-scale evaluation at Snow Crested Butte

Figure 3 presents a comparison between observations and simulations with SVS1 and SVS2 of snow height, SWE, and bulk snow density at the open site and below the forest at Snow Crested Butte. SVS2 accurately reproduced the multi-year evolution of bulk snowpack properties with error metrics that are systematically improved compared to SVS1 for the open and forested sites (Table 3). In open terrain, SVS2 simulates more accurately the snow melt dynamics and the date of snow disappearance (Fig. 3a) compared to SVS1. Such limitation of SVS1 has been identified in Leonardini et al. (2021). In the forest, SVS2 simulated a lower peak SWE than in the open due to sublimation of intercepted canopy snow and slightly underestimated observed peak SWE for the three winters (Table 3). On the other hand, SVS1 overestimated peak SWE in the forest due to the absence of snow interception and a tendency to underestimate snow melt (Leonardini et al., 2021). Finally, SVS2 better captured the increase in bulk snow density with time than SVS1 leading to improved error metrics (Table 3).

Figure 4 Simulated (with SVS1 in blue and SVS2 in green) and observed vertical snow density profiles in (a, b, e, and f) the open site and (c, d, g, and h) the forest site at two different dates during the 2018/19 winter season at Snow Crested Butte. The first row shows raw vertical profiles whereas the second row shows normalized profiles where the thickness of the layers in each observed and simulated profile is normalized by total snow height.

Figure 4 compares the simulated and observed vertical profiles of snow density for two dates of winter 2018/2019. SVS2 performed well in simulating increased snow density with snow depth, a feature that the one-layer scheme in SVS1 could not capture by design. Nonetheless, SVS2 presented an underestimation by ∼ 15 % of snow density in the top third of the snowpack in the forest compared to all the snow pit measurements. This underestimation was less present in open terrain (∼ 6 %). It may be associated with reduction of wind speed in the forest compared to the open site, which tends to reduce the snowfall density as found in Bonner et al. (2022a). Uncertainties also exist in the density of unloaded snow, which has a fixed value of 200 kg m⁻³ in SVS2/Crocus. Bouchard et al. (2024) have shown how an age-dependent density of unloaded snow can improve simulations of snow density in forested terrain. On the other hand, compared to all the snowpits over the 3 year study period, SVS2 tended to overestimate by ∼ 6 % and ∼ 10 % the density in the bottom third of the snowpack at open and forest sites, respectively.

Figure 5 (a, b) Simulated snow grain type at the (a) open site and (b) the forest site. (c, d) show the observed vertical layer grain types from the snowpits against the simulation at the open site on 16 February 2019, and 16 March 2019, respectively. (e, f) are the same as (c, d) but for the forest site. The vertical lines in (a, b) represent the times at which the comparison between the snowpits and the simulation are done (c–f). Grain type abbreviations correspond to the international classification for seasonal snow on the ground (Fierz et al., 2009).

Simulated grain type of each layer was qualitatively compared to those from snow pits (Fig. 5). No quantitative comparison was carried out due to uncertainties in the observed snow type and in the snow grain diagnostic from the model (e.g. Viallon-Galinier et al., 2020). Until mid-January, the model simulated a higher prevalence of faceted crystals and depth hoar layers in the forest relative to the open site. These results are consistent with Bouchard et al. (2022) in the humid boreal forest of Eastern Canada. In the open, the simulated snowpack presented several layers made of refrozen melt forms (MF), which were not found in measurements. A more detailed analysis showed that refrozen melt forms were simulated in the open due to several mid-winter melting events that were not simulated in the forest. Comparison with observations in February and March 2019 also showed that observed layers of faceted crystals (faceted crystals, FC, and depth hoar, DH) at the bottom of the snowpack in both the open and the forest were not fully captured by the model that simulated a succession of layers made of rounded grains (RG) and faceted crystals (FC) instead. Finally, in the forest, the model simulated more fresh snow (precipitation particles, PP, and decomposed and fragmented particles, DF) at the top of the snowpack than observed (Fig. 5e and f). This can be potentially explained by the fixed microstructure properties of unloaded snow from the canopy in SVS2 (SSA of 32 m² kg⁻¹ and sphericity of 0.95). Analogous to the treatment of density, Bouchard et al. (2024) proposed incorporating a dependency on the age of unloaded snow when computing its microstructural characteristics. This approach will be considered in a future version of SVS2.

Figure 6 Distributions of error metrics (Bias on the left and RMSE on the right) for the simulation of the vertical profile of snow density at different locations across Eastern Canada (Fig. 2) in February 2008. Error metrics are computed for the whole snowpack (a, d), the top 25 % (b, e) and the bottom 25 % (c, f). The numbers in italic at the bottom of a) and d) indicate the number of snow pits considered within each region for model evaluation (open terrain/forested terrain).

Figure 7 Same as Fig. 6 for the simulation the vertical profile of snow density for winter 2015/2016.

4.2 Distributed snowpack simulations across a latitudinal transect in Eastern Canada

4.2.1 Simulations of snowpack density

Figures 6 and 7 illustrate the performance of two SVS2/Crocus configurations (Def_LAP and Arc_LAP) in reproducing vertical snow density profiles for the winters of 2007/08 and 2015/16, respectively. The default configuration (Def) is not shown on these figures since it gave similar results as Def_LAP. Def_LAP led to a systematic underestimation of the upper snowpack densities for all locations across Quebec for winter 2007/08 (Fig. 6b). The Arctic configuration (Arc_LAP) strongly improved model results. It removed most of the negative bias in the density of the upper snowpack and decreased the median RMSE (Fig 6e) for the sub-Arctic locations (KU and QTU) but also for sites located below the treeline (SI, TQS, SC and TQN; Fig. 2). These improvements result from the modifications in the Arctic configuration that increase the compaction of surface snow due to high wind speeds and that are active over the whole simulation domain (Sect. 2.2.3). Similar improvements for sites below the treeline were found at Foret Montmorency (FMO) and Baie James (BJA) in 2015/16 (Fig 7b and d). No improvements were found at Kuujjuarapik (KUJ) in 2015/16 where both configurations underestimated the density of the upper layers. At KUJ, the snow pits were collected in a sparse forested environments where the wind speed was reduced by the presence of the trees (Eq. 17), explaining why the modifications to enhance wind packing in Arc_LAP are not as efficient as in open terrain. The Arc_LAP configuration enhanced surface snow density simulations (as evidenced by lower median RMSE values) across all sub-Arctic and Arctic sites (UMQ, CAB, RES, and WHU) during winter 2015/16 when compared to the Def_LAP experiment. However, notable simulation errors remained in the Arc_LAP results for high Arctic locations, including an overestimation of surface snow density at Cambridge Bay (CAB) and significant underestimation at the two polar desert sites (RES and WHU).

The density of upper snow layers at the sub-Arctic and Arctic sites are in agreement with previous findings of Royer et al. (2021b), Lackner et al. (2022), Woolley et al. (2024) and Woolley et al. (2025). These studies have shown that (i) the default version of SVS2/Crocus strongly underestimated the density of surface snow in the Arctic, (ii) modifications in Arc_LAP that increase the effects of wind-packing strongly improved the simulated density of upper snow layers, and (iii) the parameterization of below-canopy wind speed in SVS2 inhibits the ability of wind effect modifications to compact surface layers of the snowpack in areas of low canopy density. Transferability of parameterizations across Arctic locations presents significant challenges, as demonstrated by SVS2/Crocus simulation results. At Cambridge Bay (CAB), Arc_LAP using the Royer et al. (2021b) parameterization severely overestimated near-surface snow density, consistent with Woolley et al. (2024) at Trail Valley Creek (Northwest Territories, Canada). Polar desert environments (RES and WHU) showed the opposite trend with large density underestimation. Beyond potential ERA5 wind speed biases contributing to these errors, fundamental improvements to wind-packing parameterizations are essential. These enhancements should build on recent advances in measuring wind-induced snow transport impacts on surface snow microstructure (Walter et al., 2024) and incorporate water vapor flux representation in Crocus to capture upper layer density evolution (Domine et al., 2016b, 2018).

In the southern part of the domain (SI, TQS, SC on Fig. 6c and FMO and BJA on Fig. 7c), the Def_LAP configuration tended to underestimate basal densities. Contrasting results were obtained in this region with the Arc_LAP configuration, which improved results in 2007/08 (Fig. 6f) but led to larger errors in 2015/16 (Fig. 7f). At sub-Arctic sites, Def_LAP overestimated basal densities at KU and QTU, (Fig. 6c) and gave accurate estimation at UMQ (Fig. 7c). These sub-Arctic sites are not located in the area of Arctic vegetation (Fig. 2b) used in Arc_LAP to activate the reduction of snow compaction for basal snow (Sect. 2.2.3). As a consequence, the basal density is strongly overestimated at these sites in Arc_LAP due to overestimated wind-packing. At Cambridge Bay (CAB), located further north in the Arctic, overestimation of basal snow density is not present in Arc_LAP because basal compaction is limited by the presence of Arctic vegetation. However, despite this parameterization, the Arc_LAP configuration did not improve distributions of RMSE compared to Def_LAP at CAB (Fig 7f). This result is consistent with Woolley et al. (2024), who found no quantitative improvement in the simulation of basal snow density with their Arctic ensemble at Trail Valley Creek despite qualitative improvements in the shape of the vertical profile of snow density. Finally, in the polar desert (RES and WHU) where Arctic vegetation is not present (Fig. 2), Arc_LAP captures well the basal density at RES but underestimated it at WHU. Crocus would benefit from a representation of the water vapor flux through the snowpack to improve simulations of basal density (Domine et al., 2019) in the Arctic but also at sub-Arctic sites such as Umiujaq (Lackner et al., 2022).

The Arc_LAP configuration improved snowpack density simulations at 10 out of 13 evaluation sites in 2007/08 and 2015/16 (Fig. 6d and Fig. 7d), with improvements defined as a lower median RMSE, indicating potential for enhanced snowpack density simulation across large areas of Canada. However, further evaluation is needed in regions such the Western Cordillera and the Prairies that were not covered by the study domain considered in this paper. Evaluation at well monitored sites, e.g. ESM SnowMIP (Menard et al., 2021), are likely to reduce uncertainty associated with meteorological forcing. Remaining improvements are required in wind-packing parameterizations (Walter et al., 2024) and representation of vapor transfer in the snowpack (Domine et al., 2019).

Figure 8 Same as Fig. 6 for the simulation the vertical profile of specific surface area (in m² kg⁻¹) for winter 2015/2016. The graphics show results for the whole snowpack (a, d), the top 25 % (b, e) and the bottom 25 % (c, f).

4.2.2 Simulations of specific surface area

Figure 8 shows evaluation of two configurations of SVS2/Crocus (Def_LAP and Arc_LAP) to simulate the vertical profile of SSA for winters 2015/2016. Both configurations use the snow metamorphism scheme of Carmagnola et al. (2014) modified by Baron (2023) so that the differences in SSA are explained by differences in temperature gradients in the snowpack, resulting from differences in the snowpack density and associated thermal conductivity. Overall, the two model configurations tend to underestimate the SSA over the full snowpack (Fig. 8a). This underestimation is found in the upper layers of the snowpack (Fig. 8b) and is in agreement with Tuzet et al. (2020), who have shown a tendency of the default configuration of Crocus to underestimate near-surface SSA in alpine terrain. Arc_LAP primarily improved SSA simulations in upper layers at polar desert sites (RES and WHU), while performing similarly to Def_LAP at other locations. Woolley et al. (2024) also found improvements with the Arctic configuration for simulation of SSA in the wind slab layers at Trail Valley Creek in the Canadian Arctic.

Significant simulation errors occurred in basal snow layer SSA across all three Arctic sites (CAB, RES, and WHU; Fig. 8f). Both default configurations consistently underestimated SSA in these basal layers, which typically consist of depth hoar (Domine et al., 2018; Rutter et al., 2019; Meloche et al., 2022). The Def_LAP configuration showed latitude-dependent error increases, with median RMSE rising from 4.1 m² kg⁻¹ at CAB to 15.8 m² kg⁻¹ at WHU. While Arc_LAP improved performance by reducing median RMSE by 25 % at WHU and 55 % at RES compared to Def_LAP, with similar results at CAB, SSA underestimation in basal layers remained problematic. These findings align with Woolley et al. (2024), who reported persistent SSA underestimation in depth hoar layers for both SVS2/Crocus configurations, with only marginal RMSE improvements using the Arctic setup. These results show that the Crocus metamorphism scheme (Brun et al., 1992; Carmagnola et al., 2014) fails to simulate the evolution of depth hoar SSA. Using the parametrization of Marbouty (1980), Crocus simulates geometric growth of depth hoar, causing continuous SSA decrease to unrealistically low values. The scheme cannot distinguish between crystal geometric size of depth hoar and microscopic features that actually controls depth hoar SSA (Taillandier et al., 2007). Improvements in the simulations of the SSA of depth hoar by Crocus are required in the context of SWE retrieval and data assimilation of Ku-band backscatter in Arctic environments since depth hoar strongly scatters the radar signal (King et al., 2018; Rutter et al., 2019; Montpetit et al., 2024; Woolley et al., 2025).

Figure 9 (a) Snow melt-out date derived from the GMASI product for winter 2007/2008; (b) Difference of snow melt-out date between the Def configuration of SVS2/Crocus and GMASI; (c) Same as (b) for the Def_LAP configuration of SVS2/Crocus, (d) Same as (b) for the Arc_LAP configuration of SVS2/Crocus; (e) Same as (b) for ERA5L and (f) Same as (b) for ERA5_CRO. Results over Greenland are not shown since Greenland was masked out from the simulation domain.

Figure 10 (a) Mean snow melt-out data as a function of latitude for different snow products; (b) Mean RMSE as a function of latitude computed using GMASI as the reference and (c) forest fraction as a function of latitude. Points fully covered by water, lake and/or glaciers were removed from the analysis.

4.2.3 Simulation of snow-melt out date

Figure 9a shows the SMOD map derived from GMASI for winter 2007/2008 and Fig. 10a shows how it evolves on average with latitude. SMOD derived from GMASI ranges between 20 March-10 April (Day Of the Year, DOY, 80–100) in the southern part of the domain (south of 47° N) and the average SMOD increases continuously with latitude (Fig. 10a). Above 66° N, the average SMOD exceeds DOY 175 (24 June). Below 55° N, where the fraction of forested pixels is larger than 0.5 (Fig. 10c), the three SVS2/Crocus simulations (Fig. 9b to d) and the ERA5-Land product (ERA5L, Fig. 9e) provide similar results and overestimate the SMOD compared to GMASI. Overestimation can be explained by an underestimation of SMOD by GMASI in this region due to limitations in GMASI for detection of snow on the ground in the boreal forest as explained by Romanov (2016). Below 55° N, ERA5-Crocus (Fig. 9f) simulates a significantly lower SMOD than the SVS2/Crocus experiments and ERA5-Land (Fig. 10a). Lower SMOD is associated with the configuration of SURFEX/Crocus used in ERA5-Crocus dataset (Ramos Buarque et al., 2025) that only simulates snow in open terrain. In contrast, both SVS2 and ERA5-Land simulate snow accumulation beneath high vegetation, which tends to prolong snow cover duration in forested areas relative to adjacent open terrain under the cold climatic conditions of the boreal forest in Eastern Canada (Lundquist et al., 2013). For these reasons (limitations of GMASI in forested terrain and differences in representation of forest between the snow products), quantitative evaluation of the SMOD is only carried out above 55° N where the fraction of forested pixels is lower than 0.5 (Fig. 10c).

Table 4 Error metrics (Bias and RMSE in days) compared to GMASI for the SMOD simulations above 55° N for each SVS2/Crocus configuration and for two benchmarks products (ERA5L and ERA5_CRO) for winter 2007/2008 and winter 2015/2016.

For the SMOD evaluation, mean error metrics (bias and RMSE) compared to GMASI above 55° N were obtained for each model experiment by computing average of the SMOD values by 1° latitude bands (Fig. 10b) weighted by the number of valid grid points withing each band. Results are summarized in Table 4. Above 55° N, the default configuration of SVS2/Crocus (Def) presents a systematic underestimation of the SMOD, meaning that simulated snow cover disappears too early in this configuration (average bias of −18 d above 55° N). ERA5_CRO, which uses the same configuration of Crocus as Def, also presents a strong underestimation (average bias of −21 d above 55° N). Def and ERA5_CRO present similar spatial patterns of SMOD differences with GMASI, with a strong SMOD underestimation found on the western part of Baffin Island and on the southern part of Victoria Island (Fig. 9b and 9e). Therefore, the same Crocus configurations, implemented in two land surface schemes (SVS2 and SURFEX), driven by identical meteorological forcing and run independently, produced highly consistent results in open terrain. This supports the local-scale technical validation of Crocus within SVS2 previously conducted by the authors and highlights the model's low sensitivity to the underlying soil–vegetation scheme and associated databases in open terrain for integrated diagnostics.

The revised albedo parameterization in configuration Def_LAP led to a later snow disappearance than Def and strongly improved the performance when compared to GMASI: reduction of 50 % in the negative bias (Fig. 10a) and reduction of 38 % in the RMSE above 55° N (Fig. 10b) compared to Def. Above 55° N, Def_LAP uses values of the snow darkening coefficient, γ, larger than 150 d corresponding to an area of low deposition of LAP on snow (Fig. S5b). These values of the snow darkening coefficient are larger than the default value of 60 d used in Def leading to slower rates of decrease of snow albedo in the visible range. Gaillard et al. (2025) has shown how larger values of γ improved simulations of snow albedo in the Canadian Arctic. The results presented in Figs. 9 and 10 show how improvements in the simulation of snow albedo translate into improvements of simulated snow melt dynamics and the associated SMOD. Future enhancements to the SVS2 albedo scheme should consider how Arctic shrubs modify winter surface albedo (Belke-Brea et al., 2020) and alter solar radiation transmission within the snowpack (Domine et al., 2025).

The Arctic configuration of SVS2/Crocus (Arc_LAP) further improved the SMOD simulations. Compared to Def (Def_LAP) above 55° N, it reduced the negative bias in SMOD by 88 % (78 %) (Table 4 and Fig. 10a) and improved the RMSE by 62 % (38 %) (Table 4 and Fig. 10b). The differences compared to Def_LAP were explained by a larger peak SWE above treeline in the Arctic in Arc_LAP compared to Def_LAP as shown in Fig. 11. Larger peak SWE led to a later snow disappearance in Arc_LAP than Def_LAP, which improved estimation of SMOD compared to GMASI. The larger peak SWE was associated with lower mass loss due to blowing snow sublimation in Arc_LAP caused by denser surface snow (Fig. 6) reducing the occurrence of blowing snow events simulated by SVS2/Crocus. These results illustrate how improvements in simulation of snowpack internal properties can translate into improvements in simulation of the snowpack seasonal evolution. Results for winter 2015/2016 show similar improvements with the Arctic configuration (Table 4 and Figs. S7 and S8).

Figure 11 (Left) Maximum SWE for winter 2007/2008 simulated by experiment Arc_LAP, (Right) Relative difference of maximum SWE between experiments Def_LAP and Arc_LAP

The configuration Arc_LAP strongly improved SMOD simulations compared to two benchmarks products. RMSE was reduced by 41 % (24 %) in 2007/2008 (2015/2016) compared to ERA5L and 64 % (62 %) when compared to ERA5_CRO (Table 4). Mudryk et al. (2024) identified ERA5L and ERA5_CRO as the two top performing products for SWE climatology, variability, and trends in the Northern Hemisphere across mountain and non-mountain regions. Results presented here suggest that ERA5_CRO product could be further improved by including updated parameterizations developed in SVS2/Crocus (revised albedo scheme and new options for the Arctic snowpack). This will be easily achieved in the near future using the externalized version of Crocus that allows efficient sharing of code development between the different LSMs.

5 Conclusions

This paper presents the Soil Vegetation and Snow version 2.0 (SVS2) land surface scheme, a major upgrade to the existing SVS model (Husain et al., 2016; Alavi et al., 2016; Leonardini et al., 2021) with significant enhancements to snowpack-related processes. SVS2 integrates the externalized version of the Crocus detailed snowpack scheme (Brun et al., 1992; Vionnet et al., 2012; Lafaysse et al., 2017). Within each grid cell, SVS2 distinguishes between snow in open terrain and forested environments (if present) and simulates the snow cover and its internal properties (including snow microstructure) separately. The forested areas simulations benefit from a new canopy module in SVS2. Several innovations were also introduced in SVS2/Crocus: a revised snow albedo parameterization incorporating local climatology of LAP deposition (Gaillard et al., 2025) and new parameterizations for Arctic snowpack properties (Woolley et al., 2024).

Point-scale evaluation at Snow Crested Butte (Colorado, USA) demonstrated SVS2 capability to accurately simulate bulk snowpack properties in both forested terrain and adjacent clearings. Distributed simulations across a latitudinal transect in Eastern Canada showed that Arctic modifications of Royer et al. (2021b) and Woolley et al. (2024) substantially improved near-surface density simulations. The revised albedo scheme, combined with the Arctic configuration, reduced errors in snow melt-out date predictions by more than 50 % for winters 2007/2008 and 2015/2016 compared to the default model version, outperforming two reference products ERA5-Land (Muñoz-Sabater et al., 2021) and ERA5-Crocus (Ramos Buarque et al., 2025). However, simulations of basal layer density in the Arctic showed no quantitative improvement, suggesting the need to incorporate vapor transfer processes (Domine et al., 2019). Low model performances were also found for the specific surface area of basal snow layers.

SVS2/Crocus has demonstrated versatility across diverse snow environments (Vionnet et al., 2022; Leroux et al., 2023; Woolley et al., 2024; Gaillard et al., 2025; Meloche et al., 2024; Woolley et al., 2025). In preparation for the Terrestrial Snow Mass Mission (Derksen et al., 2021), current development efforts focus on coupling SVS2/Crocus with the Snow Microwave Radiative Transfer Model (SMRT, Picard et al., 2018; Meloche et al., 2024), the Multiple Snow Data Assimilation system (Alonso-González et al., 2022) and the Canadian Land Data Assimilation System (Carrera et al., 2015). Future developments will address several key areas: revision of the snow metamorphism scheme to improve simulations of SSA and backscatter when coupled to SMRT; implementation of water vapor flux parameterization in collaboration with the IVORI project (Brondex et al., 2023); and enhancement of the wind-packing parameterization (Walter et al., 2024). These modifications aim to establish a unified model configuration applicable across broad geographic regions, enhancing the scalability of Crocus for large-scale applications within SVS2 and other land surface models. Additionally, optimization of the maximum snow layer count will facilitate integration into operational land surface and river prediction systems (e.g., Durnford et al., 2021).