Impact of intercepted and sub-canopy snow microstructure on snowpack response to rain-on-snow events under a boreal canopy

Rain-on-snow events can cause severe flooding in snow-dominated regions. These are expected to become more frequent in the future as climate change shifts the precipitation from snowfall to rainfall. However, little is known about how winter rainfall interacts with an evergreen canopy and affects the underlying snowpack. In this study, we document 5 years of rain-on-snow events and snowpack observations at two boreal forested sites of eastern Canada. Our observations show that rain-on-snow events over a boreal canopy lead to the formation of melt–freeze layers as rainwater refreezes at the surface of the sub-canopy snowpack. They also generate frozen percolation channels, suggesting that preferential flow is favoured in the sub-canopy snowpack during rain-on-snow events. We then used the multi-layer snow model SNOWPACK to simulate the sub-canopy snowpack at both sites. Although SNOWPACK performs reasonably well in reproducing snow height (RMSE = 17.3 cm), snow surface temperature (RMSE = 1.0 °C), and density profiles (agreement score = 0.79), its performance declines when it comes to simulating snowpack stratigraphy, as it fails to reproduce many of the observed melt–freeze layers. To correct for this, we implemented a densification function of the intercepted snow in the canopy module of SNOWPACK. This new feature allows the model to reproduce 33 % more of the observed melt–freeze layers that are induced by rain-on-snow events. This new model development also delays and reduces the snowpack runoff. In fact, it triggers the unloading of dense snow layers with small rounded grains, which in turn produces fine-over-coarse transitions that limit percolation and favour refreezing. Our results suggest that the boreal vegetation modulates the sub-canopy snowpack structure and runoff from rain-on-snow events. Overall, this study highlights the need for canopy snow property measurements to improve hydrological models in forested snow-covered regions.

1 Introduction

There are several definitions of a rain-on-snow (ROS) event (Brandt et al., 2022). Simply put, an ROS event occurs when liquid precipitation falls on an existing snowpack (McCabe et al., 2007). As simple as the definition of an ROS event may be, these can have major consequences if the right conditions are met (Wayand et al., 2015). Historically, some winter ROS events have caused severe flooding and damage in North America, Europe, and many other regions (Marks et al., 1998; Haleakala et al., 2023; Rössler et al., 2014). The risk of flooding from ROS events is controlled not only by atmospheric and soil conditions (Haleakala et al., 2023; Zaqout et al., 2023), but also by the energy balance and structure of the snowpack. A warm and thin snowpack favours a greater contribution of snowmelt to flooding (Jennings et al., 2018; Brandt et al., 2022), while a cold and thick snow cover limits the occurrence of snowmelt (Trubilowicz and Moore, 2017). In cold and stratified snowpacks, preferential flow tends to occur, causing rapid runoff during ROS events (Avanzi et al., 2016; Würzer et al., 2017). In contrast, snowpack runoff is slowed by ice layers and melt–freeze crusts within the snowpack that impede water percolation (Webb et al., 2018a; Würzer et al., 2016; Eiriksson et al., 2013). In the Northern Hemisphere, ROS events are expected to become more intense with warmer winters in response to a shift from solid to liquid precipitation, increasing the risk of flooding (IPCC, 2023; Musselman et al., 2018).

The relationship between the presence of a forest and runoff from ROS events remains unclear, as only a few studies have explored this topic (Brandt et al., 2022). Beaudry and Golding (1983) and Berris and Harr (1987) found that tall and dense coniferous forests in western North America experienced reduced runoff during ROS events when compared to open sites. In contrast, Berg et al. (1991) and Garvelmann et al. (2015) found that vegetation has a small impact on runoff following an ROS event in mature and dense mixed forests of the western United States and central Europe, respectively. None of these studies have examined the role of the sub-canopy snowpack structure during ROS events, despite its strong influence on water transport and snowpack runoff.

Vegetation heavily influences the evolution of snowpack structure in snow-dominated regions. During winter, evergreen canopies intercept a substantial fraction of the solid precipitation (Hedstrom and Pomeroy, 1998; Pomeroy et al., 2002), limiting snow accumulation underneath (Varhola et al., 2010; Parajuli et al., 2020). A thin snowpack below the canopy increases the vertical temperature gradient, favouring gradient metamorphism and grain growth (Musselman et al., 2008; Molotch et al., 2016). Confirming this, Bouchard et al. (2022) observed that the snow density under the canopy is lower and the grains are larger and have a lower specific surface area (SSA) than for the snowpack within the adjacent forest gaps. The combination of a lower density and a lower SSA leads to a greater permeability of the sub-canopy snowpack, facilitating the downward water flow compared to the gaps. Previous studies have also shown that preferential flow is more likely under the canopy due to canopy snow unloading, meltwater dripping, and accumulation of vegetation debris over the sub-canopy snow cover (Bründl et al., 1999; Teich et al., 2019). Given the influence of forest canopies on snowpack structure, one can expect that the runoff generated from ROS events would also be impacted by the presence of vegetation. However, studies are needed to confirm this.

Detailed snow models can help capture the complex dynamics of water flow through the snowpack and aid in the interpretation of field observations. SNOWPACK (Lehning et al., 2002a) is a multi-layer, one-dimensional snow model that can solve Richards' equation to describe water percolation from preferential flow within the snow cover (Wever et al., 2014, 2016). SNOWPACK uses a big-leaf approach to represent the canopy (Lehning et al., 2006). A two-layer canopy module with thermal inertia has been implemented, which allows better representation of the diurnal variation in canopy temperature and heat exchanges with the sub-canopy snowpack (Gouttevin et al., 2015). SNOWPACK has been evaluated in forested environments with regard to snowpack energy balance (Gouttevin et al., 2015; Todt et al., 2018), snow accumulation (Rutter et al., 2009; Gouttevin et al., 2015), and snow-related climate feedbacks (Krinner et al., 2018). However, only Rasmus et al. (2007) and Kontu et al. (2017) used detailed snow profile observations to validate SNOWPACK in the boreal forest of Finland. In both studies, the authors found a reasonable agreement between the model and the observations, but their model validation was limited to open forest clearings, without consideration for snow under the trees. Further evaluation of SNOWPACK to simulate the sub-canopy snow cover in northern regions is therefore in order.

SNOWPACK, as any other snow model, tends to underperform in forests when compared to open areas (Rutter et al., 2009). This is partly due to a lack of calibration and validation snow data in forests and to interception parameterizations derived from a few observational studies that are difficult to generalize to other climates (Lundquist et al., 2021). In SNOWPACK, the interception scheme is based on experiments by Schmidt and Gluns (1991) and Hedstrom and Pomeroy (1998) conducted in a continental climate of the western United States and Canada. In this parameterization, the intercepted fraction of the precipitation is stored in the canopy until the temperature-dependent maximum capacity is reached (Pomeroy et al., 1998). The intercepted snow is then unloaded when the canopy storage exceeds its maximum capacity. The unloaded snow is considered “fresh snowfall”, with properties defined as a function of meteorological variables such as air and snow surface temperature, relative humidity, and wind speed at the time of the unloading event (Lehning et al., 2002b). However, snow can remain in the canopy for up to several weeks (MacDonald, 2010; Lumbrazo et al., 2022) and undergo metamorphism before being unloaded. Since these parameterizations do not account for snow metamorphism, this can lead to inaccurate simulations of the physical properties and stratigraphy of the sub-canopy snowpack after an unloading event. This, in turn, could alter the simulated downward water flow from an ROS event. Moreover, the high spatial heterogeneity of the unloading mechanisms makes it even more difficult to accurately simulate the sub-canopy snowpack using a one-dimensional snow model (Vincent et al., 2018).

In this study, we address ROS events in the boreal forest through observations and modeling. Specifically, we first use field observations to document the ROS-induced alterations of sub-canopy snowpack structure, with a focus on melt–freeze snow layers and preferential flow channels. This objective builds on the work of Bouchard et al. (2022), who found that the sub-canopy snowpack is highly heterogeneous and has high permeability, expected to facilitate downward water flow. These field observations provide the unique opportunity to evaluate the one-dimensional, multi-layer, microstructure-resolving model SNOWPACK in the boreal forest of eastern Canada. Motivated by the hypothesis that unrealistic properties of unloaded snow hamper the performance of SNOWPACK in simulating how sub-canopy snowpacks react to ROS events, we suggest an alternative representation of canopy snow properties and assess its impact on SNOWPACK simulations of melt–freeze layers and ROS-induced runoff. Overall, this work improves both our understanding of ROS events in the boreal forest and our ability to capture their impacts on snowpack structure and runoff thanks to physically based modeling. This is especially important as ROS events are expected to become more frequent in the future.

2 Observational data

2.1 Study sites

In this study, we compare two sites located in the boreal forest of eastern Canada where meteorological variables are continuously measured above the canopy. The Montmorency Forest (MF) is the main site, as it is more thoroughly instrumented and readily accessible in winter. The Bernard River Valley (BRV) is located some 700 km northeast of MF and requires a full day's drive and an hour ski-in to reach in winter. This site, which covers a different bioclimatic area than MF and where less data have been collected, is therefore used primarily to validate conclusions drawn from the main site.

2.1.1 Montmorency Forest

MF (47°17^′18^′′ N, 71°10^′05^′′ W) is located in the province of Quebec, Canada, at the southern edge of the boreal forest (Fig. 1a). The site receives an average (1991–2020) total annual precipitation of 1504 mm, of which about 40 % falls as snow (station no. 7042395 Foret Montmorency, https://climate.weather.gc.ca/climate_normals/index_e.html, last access: 9 October 2023). The mean 2 m air temperature in December, January, and February (DJF) is −13.2 °C.

Figure 1 Study sites. (a) Location in the province of Quebec (QC), Canada. (b, c) Hemispherical image of sub-canopy locations at both sites along with an artistic depiction of the dominant tree species at each site (balsam fir – MF; black spruce – BRV). (d–f) Photos taken from the meteorological towers showing the forest canopy. (e–g) Sub-canopy monitoring stations showing the thermistors array, snow height sensor, and air-temperature–relative-humidity probe at both sites.

The field campaign was conducted on a 12° northeast-facing slope at an elevation of ≈850 m above mean sea level (a.m.s.l.) within MF. The site is a juvenile balsam fir (Abies balsamea) and white birch (Betula papyrifera) stand that naturally regrew after harvesting operations in 1993–1994 (Guillemette et al., 2005). The canopy height is between 7 and 12 m, with an average height of 9.2 m when estimated in the summer of 2019 using a clinometer. The soil is a sandy loam overlaid by a 7 cm layer of organic material. The average leaf area index (LAI) of sub-canopy locations at MF is 4.8±1.6 m² m⁻². The LAI was derived from the analysis of 26 hemispherical canopy images using the method from van Gardingen et al. (1999), designed for highly clumped canopies such as dense coniferous stands. In this approach, the image is separated into 90×1° zenith rings and the LAI is calculated from the non-linear estimation of the log-averaged gap fraction of each ring. In the summer of 2015, a 15 m flux tower was erected on site to measure hydrometeorological variables (Isabelle et al., 2018). In the fall of 2020, the tower and the instruments were raised to a height of 20 m. Air temperature (T_a) and relative humidity (RH) were both measured using an HC2S3 probe (Rotronic). Incoming and outgoing shortwave (SWR_↓ and SWR_↑) and longwave radiation (LWR_↓ and LWR_↑) were monitored with a CNR4 radiometer (Kipp & Zonen) inclined so radiation was measured perpendicular to the slope. Wind speed (WS) was measured using a two-dimensional anemometer (R.M. Young – model 05103). Table A1 shows the measurement height of each meteorological variable at MF before and after raising the tower in October 2020. All variables are sampled every minute and averaged to the half-hour.

Finally, a double fence automatic reference (DFAR) for liquid and solid precipitation (P) is located in an open clearing about 4 km northeast of the study site at ≈670 m a.m.s.l. (Pierre et al., 2019). We assume that P measured by the DFAR is representative of the precipitation received at the experimental site despite the distance and the elevation difference.

2.1.2 Bernard River Valley

The BRV site (50°54^′36^′′ N, 63°22^′48^′′ W) receives 1077 mm of precipitation annually and the DJF temperature is −12.5 °C, according to the 1991–2020 climate averages from the nearest federal weather station, located about 200 km southwest from the site (station no. 7047914 Sept-Iles, https://climate.weather.gc.ca/climate_normals/index_e.html, last access: 9 November 2023). The valley has a mean elevation of ≈250 m a.m.s.l. and is surrounded by plateaus 200 m higher in elevation. BRV is a site dominated by black spruce (Picea mariana), with trees between 12 and 18 m high, with an average height of 15 m, when estimated in the spring of 2021, also with a clinometer. The soil is composed of silty loam mineral soil under an organic layer of 17 cm. Compared to MF, the forest at BRV is rather sparse (Fig. 1c and f). We estimated an LAI of 1.6±0.7 m² m⁻² from 55 LAI-2200C (Li-COR) measurements, also taken in the spring of 2021. A 25 m flux tower is used to measure the same meteorological variables as at the MF site. Air temperature and relative humidity were measured using an HMP device, all-wave radiation was measured by a CNR4, and wind speed was obtained from two-dimensional anemometer measurements. The same sampling frequency and averaging time step as at MF are used, and the height of the instrument is presented in Appendix A.

2.2 Snow observations

Snow observations cover the period of 2018 to 2023 at MF (five winters) and 2019 to 2023 at BRV (four winters). Sub-canopy monitoring stations were installed approximately 25 m from the meteorological tower at each site prior to the first field campaign. Continuous automated measurements were complemented by recurrent snow pit observations. These were conducted at sub-canopy sites with similar characteristics to the monitoring station and within a radius of 100 m of the tower at both sites.

2.2.1 Snow monitoring stations

At the MF site, snow height was monitored with an SR50 ultrasonic sensor (Campbell Scientific) in winter 2018–2019 and with a Judd communication snow height sensor from 2019 to 2023. Snow surface temperature (T_surf) was monitored with an SI-111 infrared radiometer (Apogee Instruments) in 2018–2019 and from 2020 to 2023. Snow height and T_surf measurements were recorded at an hourly time step on a CR10X data logger (Campbell Scientific) from mid-October to mid-June. During winters of 2020–2021, 2021–2022, and 2022–2023, we also monitored snowpack temperature every 15 cm from the ground level with Pt-1000 thermistors (Schneider Electric) and T_a at 2.5 m height using an HMP probe (Fig. 1e and g). A more detailed description of the automatic stations can be found in Bouchard et al. (2023a).

At the BRV sites, snow height was measured using an SR50 ultrasonic sensor from 2019 to 2021 and then with a Judd communication snow height sensor from 2021 to 2023. We also monitored snowpack temperature every 15 cm using Pt-1000 thermistors in a setup similar to that described for MF. T_surf was not measured at BRV.

At both the MF and BRV sites, we installed time-lapse cameras on trees at a height of 2.5 m and facing the canopy. This setup allowed us to detect the presence of snow in the canopy and to identify the phase of the precipitation through visual interpretation of the time-lapse images.

2.2.2 Snow pit measurements

A total of 48 snow pits were dug under the canopy during the 2018–2023 study period (42 at MF and 6 at BRV). Of this number, 26 snow pits were dug at MF during the winter of 2018–2019 (Bouchard et al., 2022). All snow pits included a visual assessment of the snowpack stratigraphy, with detailed identification of ice and melt–freeze layers within the snowpack. Grain types were identified using a magnifying glass and a millimetric grid by the same observer throughout the study period. Snow pits also included a measurement of the vertical profile of snow density, except for the last snow pit of winter 2022–2023 at MF. Snow density was measured every 3 to 5 cm with a 100 cm³ density cutter box with a ±10 % accuracy (Conger and McClung, 2009).

2.3 Soil observations

Soil profile characterizations were performed in the fall of 2020 (BRV) and the summer of 2021 (MF). We measured a vertical profile of temperature using a Pt-1000 thermistor (Senseca), soil thermal conductivity using a TP02 heated needle probe (Hukselflux), and soil volumetric water content and density by gravimetric analysis, along with a soil texture identification through the profile. Soil properties were measured up to depths of 100 cm at MF and 82 cm at BRV. More details on the observed soil characterization are provided in the Supplement (Figs. S1 and S2).

2.4 Rain-on-snow events

In this study, we define an ROS event as at least 3 mm of rain cumulated over 12 h or longer while a minimum of 3 cm of snow on the ground is observed. It is one of the many definitions of an ROS event in the literature (Brandt et al., 2022), which we chose for its simplicity in the absence of runoff measurements. Rainfall amount and duration were taken from the DFAR at MF and from ERA5-Land reanalysis at BRV (see Sect. 3.4), whereas precipitation phase and snow height were respectively obtained from time-lapse images and ultrasonic measurements at both sites. Note that we focus on ROS events that occur between November and March exclusively.

3 Numerical modeling

We used the one-dimensional physically based snow cover model SNOWPACK (version 3.6.0) to simulate the sub-canopy snowpack from 2018 to 2023 in the Montmorency Forest and from 2019 to 2023 in the Bernard River Valley. The SNOWPACK model was selected for several reasons. First, it accounts for the processes that drive snow metamorphism and provides a complete and detailed representation of snow microstructure, thermal profile, snow settlement, and mass and energy balance (Bartelt and Lehning, 2002; Lehning et al., 2002a, b). Second, it solves Richards' equation to represent matrix and preferential flow in the snowpack (Wever et al., 2016). Third, SNOWPACK includes a two-layer canopy module that simulates the exchanges of mass and energy between the vegetation and the underlying snowpack (Gouttevin et al., 2015). Our detailed snow dataset provides an excellent opportunity to validate and further improve this model for forest applications. Simulations were run at a 15 min calculation time step. Contents of SNOWPACK initialization files are presented in Appendix B.

3.1 Water transport scheme

The water transport scheme in SNOWPACK includes two options: a simple bucket model approach and a solution of Richards' equation (Wever et al., 2014). The latter approach allows for better simulation of snowpack runoff volume and timing, particularly in the early melt season and on sub-daily timescales (Wever et al., 2014, 2015). Wever et al. (2016) further implemented a description of preferential flow using a dual-domain approach, where water is transferred from the preferential to the matrix flow domain when it encounters a layer transition, to simulate fine-over-coarse snow layers (Katsushima et al., 2013; Avanzi et al., 2016). In the model, this is conceptually characterized by a threshold of liquid water content that marks the saturation of the preferential flow domain (θ_TH; 0–1). In practice, the use of θ_TH controls the movement of water from preferential to matrix flow. This allows for water spreading at microstructure transitions, which leads to the formation of ice or melt–freeze layers, as phase change is only possible in the matrix flow domain. An ice reservoir parameterization was further developed by Quéno et al. (2020) to better capture the formation of continuous ice layers from discontinuous and growing ice lenses. This development improved the formation of ice and melt–freeze layers and reduced the number of simulated ice layers that were not observed.

We tuned θ_TH based on SNOWPACK simulations at an open site, i.e., without vegetation effects, where we examined snowpack wetting after an ROS event. The simulations were compared with snow pit observations performed at the open site, where thin melt–freeze and ice layers resulting from ROS events were identified (unpublished data). These were better reproduced by simulations with a θ_TH greater than 0.35. Therefore, we set θ_TH to 0.35 in our modeling setup. Moreover, the snow water equivalent appeared to be almost unsensitive to changes in θ_TH.

3.2 SNOWPACK canopy module

The canopy module of SNOWPACK relies on a two-layer canopy scheme as detailed by Gouttevin et al. (2015). The mass balance in the canopy depends on snow interception, canopy evaporation and sublimation, and unloading. The interception rate is calculated from the canopy storage saturation, which varies with intercepted snow density (ρ_s,int) and LAI. The parameterization of ρ_s,int is the same as for new snow (Lehning et al., 2002a). Evaporation and sublimation are calculated as part of the two-layer energy balance and added to the mass balance at the end of each time step. A complete description of the energy balance parametrization can be found in Sect. 2.4 and 2.5 of Gouttevin et al. (2015). Unloading occurs when canopy storage exceeds the maximum storage capacity of the canopy. At the onset of an ROS event, snow in the canopy unloads first and contributes to the formation of a new snow layer on top of the snowpack with fresh snow properties based on the weather conditions at the time of unloading. A more detailed description of interception and unloading parameterization can be found in Appendix C.

3.2.1 Intercepted snow densification

We expect SNOWPACK simulations under the canopy to be sensitive to the properties of the unloaded snow. In its current form, SNOWPACK does not simulate the evolution of snow properties in the canopy. Therefore, we implemented a canopy snow parameterization that accounts for the density and microstructure of the intercepted snow (Fig. 2). From now on, we will refer to this version of the canopy module as “ISD” for “intercepted snow densification”, while the original version of the canopy module will be referred to as the “initial module” or “IM”.

In ISD, the density of the intercepted snow (ρ_s,int; in kg m⁻³) can be estimated as a function of its age (a_s,int; in days) based on Eq. (16) from Koch et al. (2019), which was originally developed to estimate snow height along with snow water equivalent measurements from an alpine site in Switzerland:

$$\begin{array}{}\text{(1)}& {\mathit{\rho }}_{\mathrm{s},\mathrm{int}}={\mathit{\rho }}_{\mathrm{fr}}+({\mathit{\rho }}_{\max }-{\mathit{\rho }}_{\mathrm{fr}})\left(\mathrm{1}-e^{-\frac{a_{\mathrm{s},\mathrm{int}}}{\mathit{\tau }}}\right),\end{array}$$

where ρ_fr and ρ_max are the fresh snow density and the maximum density for intercepted snow (in kg m⁻³), respectively, and τ is the shape parameter of the exponential function. We chose to use the model from Koch et al. (2019) because it could be calibrated using the field data we have available.

Using snow pit density measurements, we determined ρ_fr as an average of 11 density measurements of precipitation particle (PP) layers taken during the study period at both sites, which resulted in a value of 80.3 kg m⁻³. For ρ_max, we used the average of 69 measurements of rounded-grain (RG) snow layers during the same period, which led to a value of 280.5 kg m⁻³. τ was defined such that ρ_s,int is within a 1 % deviation of ρ_max after $a_{\mathrm{s},\mathrm{int}}=\mathrm{30}$ d, as in Koch et al. (2019). If new snow is intercepted during a time step, a_s,int decreases proportionally to the amount of new snow added to the canopy, as follows:

$$\begin{array}{}\text{(2)}& a_{\mathrm{s},\mathrm{int}}=a_{\mathrm{s},\mathrm{int}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{\Delta }I\mathrm{\Delta }t}{(I+\mathrm{\Delta }I\mathrm{\Delta }t)}}\right),\end{array}$$

where I is the interception storage (in mm), ΔI is the interception rate in mm d⁻¹, and Δt is the 15 min computational time step (0.010417 d). As long as ρ_s,int is less than a threshold density (ρ_th), the snow stored in the canopy is considered “fresh snow”. This means that the microstructure parameters (dendricity (dd), sphericity (sp), grain diameter (d_g), and bond size (d_b)) are those of fresh snow according to the parameterization of Lehning et al. (2002b). Similar to before, we empirically set the ρ_th value to 152 kg m⁻³, corresponding to the average of 26 density measurements in layers of decomposed and fragmented (DF) precipitation particles taken during the study period at both sites. This threshold value is reached after 3.12 d according to Eq. (1) if no additional snow is intercepted in the meantime. When ρ_s,int exceeds the threshold, snow stored in the canopy is considered as RG with dd = 0, sp = 1, d_g=0.2 mm, and $d_{\mathrm{b}}=d_{\mathrm{g}}/\mathrm{3}=\mathrm{0.07}$ mm. Parameters dd, sp, and d_b are all default values for RG in SNOWPACK (Lehning et al., 2002b). Grain diameter d_g was set to 0.2 mm, which corresponds to a specific surface area (SSA) of 32.7 m² kg⁻¹ for an optical equivalent grain size (Grenfell and Warren, 1999):

$$\begin{array}{}\text{(3)}& \text{SSA}={\displaystyle \frac{\mathrm{6}}{{\mathit{\rho }}_{\mathrm{ice}}{d}_{\mathrm{g}}}},\end{array}$$

where ρ_ice is the density of ice (917 kg m⁻³). These values of d_g and SSA are within the range of multiple sub-canopy measurements of RG from previous studies (Molotch et al., 2016; Bouchard et al., 2022, 2023a).

Figure 2 Schematic description of the initial canopy module (IM) and the intercepted snow densification (ISD) with respect to canopy snow properties. With IM, the density of unloaded snow (ρ_u) is equivalent to the density of new snow (ρ_new) and is a function of air temperature (T_a), relative humidity (RH), wind speed (WS), and snow surface temperature (T_surf). ISD tracks the age of the intercepted snow (a_s,int). When new snow is intercepted, a_s,int is adjusted as a function of interception storage (I) and rate (ΔI) before getting older by one time step (Δt). Then, the density of the intercepted snow (ρ_s,int) is computed based on a_s,int. Snow unloads as precipitation particle (PP) snow with the corresponding dendricity (dd), sphericity (sp), grain diameter (d_g), and bond diameter (d_b) when ρ_s,int is below a threshold density (ρ_th). If ρ_s,int is larger than ρ_th, snow unloads as rounded grains (RG) with the corresponding dd, sp, d_g, and d_b. In both cases, ${\mathit{\rho }}_{\mathrm{u}}={\mathit{\rho }}_{\mathrm{s},\mathrm{int}}$. Symbols of PP and RG are taken from the International Classification for Seasonal Snow on Ground (Fierz et al., 2009).

3.2.2 Canopy module parameters

Table 1 shows the canopy parameters that we used for simulations at both sites. Canopy height (z_can; m) and LAI were based on field measurements described in Sect. 2.1. The stand basal area (B; m² m⁻²), which refers to the fraction of total surface area occupied by trunks, was taken from Hadiwijaya et al. (2020) at the MF site. We used half of this value at the BRV site due to the sparser canopy. Since our focus is on snow accumulating below the canopy, we set the direct throughfall fraction (c_f) to 0 in our simulations at both sites. For the interception capacity factor, parameter i_max, we applied the value suggested by Schmidt and Gluns (1991) for spruce at both sites. We used the wet, dry, and snow-covered canopy albedo (α_wet, rain, α_dry, and α_wet, snow) and the LAI fraction for the needle layer (f_LAI) from Gouttevin et al. (2015), who showed that these parameters were easily transferable between sites. A full description of all the canopy parameters can be found in Gouttevin et al. (2015).

Table 1 User-provided model parameters for the canopy module.

3.3 Soil parameterization

In SNOWPACK, snow and soil layers form one continuous column. Snowpack runoff is defined as the liquid water transferred from the lowermost snow layer to the uppermost soil layer and is computed by applying Darcy's law at the soil–snow interface (see Eq. 15 from Wever et al., 2014). The water content in the soil is further calculated based on the van Genuchten model (see Eq. 3 in Wever et al., 2014).

The soil was parameterized from in situ measurements (see Sect. 2.3). The temperature at the bottom of the lowest soil layer, at a depth of 200 and 182 cm from the surface for MF and BRV simulations, respectively, was defined as a Dirichlet boundary condition, where the specified temperature corresponds to the deepest measurement. Each soil layer was set with a thickness of 1 cm, and the properties of lowest observed soil layer were replicated down to the lower boundary at both sites in the simulation setup. The initial soil parameters provided to SNOWPACK are shown in the Supplement (Tables S1 and S2).

After initializing the soil, we performed a 10-year and 12-year spinup at MF and BRV, respectively, to stabilize the water content and the thermal regime of the soil (Rodell et al., 2005). In the spinup simulations and in the absence of longer measurement time series, we looped winter conditions twice from 2018 to 2023 at MF and three times from 2019 to 2023 at BRV to obtain a state of equilibrium for the soil.

3.4 Forcing data

In our study, SNOWPACK is driven in offline mode at the point scale (one-dimensional) using local meteorological observations of T_a, RH, WS, SWR_↓, LWR_↓, P, and precipitation-phase (P_ph) data and run at a 30 min time step. In the Bernard River Valley, we used the hourly precipitation from the ERA5-land reanalysis (Muñoz-Sabater et al., 2021), which we equally divided into 30 min time steps. The ERA5-Land reanalysis shows a good agreement with a Hydro-Québec gauge located 20 km east of the study site for the annual cumulative precipitation (Fig. S3). The Hydro-Québec gauge was not used as forcing data because it only measured precipitation at the millimeter resolution. Since the phase of the precipitation is used to estimate the canopy interception capacity (see Appendix C), a linear transition with dual temperature thresholds at 0 and 2 °C was first applied to define the phase of precipitation at any time step. Then, this phase was validated using time-lapse images. Note that in the model, fluxes were calculated perpendicular to the slope according to the radiation measurements above the canopy.

3.5 Model evaluation

We first evaluated the model performance in simulating snow height using the root-mean-square error (RMSE; cm) and the bias relative to the observations (p_bias; %). We also evaluated the difference in snow disappearance date (ΔSDD; days). We further calculated the RMSE of T_surf (in °C) at MF in 2018–2019 and from 2020 to 2023. Validating simulated T_surf provides an estimate of how well the overall energy balance of the snowpack is captured. Then, we compared the simulated and observed vertical profiles of density and grain type and evaluated the agreement between both on a score from 0 to 1 (agreement score) based on the method described by Lehning et al. (2001). Finally, we visually compared the simulated stratigraphy with field observations of melt–freeze and ice layers as a semi-quantitative analysis of the number of correctly simulated melt–freeze layers.

Given that ρ_fr, ρ_max, ρ_th, and d_g values can be subject to site dependencies, we evaluated the sensitivity of the ISD module to these canopy snow parameters. Table 2 shows the low and high values used in the sensitivity analysis, as well as the values assigned for the baseline analysis. For the sensitivity analysis, we manually changed the values of one parameter at a time to isolate its effect on the simulation results. We performed the sensitivity analysis exclusively for winter 2018–2019 at MF, as this is the year with the most snow pit observations. We first objectively assessed the melt–freeze layer formation, as described in Sect. 3.2 from Quéno et al. (2020). Briefly, this method evaluates whether the observed layer is well simulated (“hit”) and whether the layer thickness is adequately reproduced. The simulated and the observed layer height must agree within a margin of 20 % of the observed total snow height to be considered as a “hit”. The observed melt–freeze layer thickness is well reproduced when the simulated layer thickness is half to twice that of the observed layer. We also evaluated the sensitivity of SNOWPACK to ROS-induced runoff, although no field data were available to validate this process.

Table 2 Values of canopy snow parameters used for the sensitivity analysis.

4 Results

The Results section is divided into three parts. In the first part, we present the ROS events that took place at the MF and BRV sites, along with the snowpack observations at both sites. In the second part, we describe the evaluation of the SNOWPACK model using the IM version of the canopy module. In the third part, we present the effects of implementing the ISD parameterization into SNOWPACK on snow unloading, snowpack structure, and snowpack runoff during ROS events, along with a case study of a specific ROS event and the sensitivity analysis of the ISD function parameters.

4.1 Climatic conditions

4.1.1 Snow height and ROS events

At both sites, snow began to accumulate between late October and mid-November and started to decline between late March and early April (Fig. 3). In general, the snowpack disappeared in late May at MF and in early May at BRV, for a snow cover duration on average 15 d longer at the MF site (197 versus 182 d). Also, more snow tended to accumulate at the MF site than at the BRV site. However, winter 2020–2021 at MF was by far the year with the lowest snow accumulation among all site years. This winter was an exceptionally warm and low-snow year at MF and is analyzed in detail in Bouchard et al. (2023a).

Figure 3 Observed snow heights from October 2018 to July 2023 at Montmorency Forest (MF; a) and from November 2019 to July 2023 at Bernard River Valley (BRV; b). Vertical blue bars show the rain-on-snow (ROS) events recorded from November to March at both sites (47 in total), with darker blue indicating greater rainfall accumulation per event (see colour code). The precipitation phase was determined using time-lapse cameras. Vertical red bars indicate the dates of snow pits (48 in total).

We observed 27 ROS events at MF and 20 at BRV from November to March throughout the study period (Figs. 3 and 4). A total of 22 ROS events were less than 10 mm, and only 7 of them exceeded 30 mm. We recorded one event of more than 100 mm at each site, on 24 December 2020 at MF (106 mm) and on 1 December 2020 at BRV (113 mm). In total, 31 ROS events occurred from October to December and 16 from January to March. Few events were observed when the air temperature was below freezing, none of them being larger than 25 mm. In general, ROS events occurred when the air temperature was between 0 and 5 °C, although a few large events took place when the temperature was higher than 5 °C. Long ROS events did not lead necessarily to important rainfall accumulation, as shown by the large dots in the lower part of Fig. 4. In contrast, events of more than 30 mm all lasted longer than 12 h.

Figure 4 Total liquid water accumulation, mean air temperature, and duration of each ROS event observed for the November to March period from 2018 to 2023 in the Montmorency Forest (MF) and from 2019 to 2023 in the Bernard River Valley (BRV).

4.1.2 Snowpack observations

During the study period, 18 melt–freeze layers were identified at MF and 14 at BRV (Table 3). Melt–freeze layers formed at the surface by diurnal melt–freeze cycles are excluded from the list, as are layers formed by ROS events from April to June. Also, since the observations were made under the forest cover, it is unlikely that shortwave radiative melting followed by nocturnal longwave radiative cooling led to a melt–freeze layer (Höller, 2001; Malle et al., 2019). Therefore, it is assumed that the melt–freeze layers listed in Table 3 were all formed at the snow surface and originated from an ROS event between November and March.

Table 3 All observed melt–freeze layers at Montmorency Forest (MF; 18) and at Bernard River Valley (BRV; 14) that are assumed to have formed on the snow surface after a rain-on-snow (ROS) event. Details include the formation date^∗, the average height of the bottom of the layer, its average thickness, and the number of times each layer was observed during snow pit experiments. The cross indicates a basal melt–freeze layer.

^∗ The formation date is assumed to be the date of the first rain-on-snow event that occurred after the date of the snow pit acquisition during which the layer beneath the melt–freeze layer was observed for the first time.

At both sites, melt–freeze layers were always observed at the base of the snowpack, except for winter 2022–2023. Not surprisingly, most layers were formed from October to December, as ROS events were more frequent during that period. Therefore, more melt–freeze layers were observed in the bottom 50 cm of snow. Note that the average thickness of the melt–freeze layers was highly variable, ranging from 1 to 9 cm.

Several small frozen percolation channels between two layers (BRV22–23c and BRV22–23d from Table 3) are shown in Fig. 5a and b, indicating that preferential flow occurred through multiple small channels within a short distance. In Fig. 5c, three small percolation channels were observed above a melt–freeze layer (MF21–22b from Table 3), which was itself above one larger percolation channel. This suggests that less preferential flow paths form deeper in the snowpack where snow grains are coarser, consistent with the observations of Katsushima et al. (2013). Pockets of snow with liquid water were also observed at various depths in the snowpack (Fig. 5d). The absence of a continuous vertical wetting between these pockets suggests that water was transported downward by preferential flow. Also, since there was no discontinuity between the snow identified as “water pockets” and the surrounding snow, we believe that these were pockets of liquid water rather than pieces of unloaded snow. Figure 5e–f show ice clumps on the surface of the snowpack during snowmelt. These ice structures, which show a clear vertical shape, appear to be residual preferential channels melting at slower rates than the surrounding snow, as also reported by Teich et al. (2019).

Figure 5 (a) Melt–freeze layers and percolation channels (30 March 2023 – BRV). (b) Closer view of the leftmost percolation channel observed in (a). (c) Preferential flow paths are less numerous but larger in diameter deeper in the snowpack where the snow grains are coarser (10 March 2022 – MF). (d) Pockets of liquid water within the snowpack (30 March 2022 – BRV). (e–f) Residual percolation channels remaining on top of the sub-canopy snowpack during the snowmelt period (28 April 2022 – MF).

Figure 6 shows observations of air temperature, snow temperature, and temperature at the soil–snow interface for an ROS event at each site. In the first event (Fig. 6a–c), air temperature rose to 0 °C on the night of 15–16 February 2023 at MF, before falling to −23 °C over the next 2 d. A total of 11.5 mm of rain was measured over 11 h when the air temperature was highest. The top layers of the snowpack reached 0 °C quickly after the start of the ROS event. The snow underneath gradually warmed from above but never reached 0 °C, except for the temperature probe at 15 cm height, which reached 0 °C some 24 h later. This vertically inhomogeneous warming pattern suggests preferential flow. A constant soil–snow temperature around −0.3 °C during the ROS event suggests that percolation did not reach the ground.

In the second event (Fig. 6d–f), a temperature rise above 0 °C at BRV on 6 and 7 March 2023 triggered a phase change in the precipitation, in which 5.3 mm of rain was recorded during this event that lasted for almost 40 h. This ROS event unfolded in two main steps. First, rain and warm air warmed the top snow layers to 0 °C, from the surface down to a snow height of 45 cm. Then, additional rain led to a rapid increase in the soil–snow interface temperature to 0 °C before snow layers above. This supports the idea that water reached the ground through preferential flow. Note that given the important water percolation in the snowpack from that event, it is likely that ERA5 underestimated rainfall accumulation.

Figure 6 Evolution of the 2.5 m high air temperature (a, d), snowpack temperature profile (b, e), and soil–snow interface temperature (c, f) during and after the ROS event of 16 February 2023 at Montmorency Forest (MF; a–c) and during the ROS event of 6 March 2023 at Bernard River Valley (BRV; d–f). Liquid precipitation and solid precipitation are shown in panels (b) and (e), and gray-red areas in the temperature profile indicate a temperature of 0 °C. Snowpack temperature is vertically interpolated between measurement heights (every 15 cm), and the observed temperature was forced to a maximum of 0 °C. Snow height is shown by the black line in panel (e) since snow surface temperature was not measured at BRV.

4.2 Evaluation of SNOWPACK

SNOWPACK is hereby evaluated with the initial interception module (IM). SNOWPACK with IM overestimates the snow height under the canopy at MF with a RMSE of 17.3 cm and a p_bias of 19.5 % (Fig. 7). Also, the simulated snowpack melts out later than the observations (ΔSDD = 5.5 d). The model performance is lower at BRV than at MF, which could be partly explained by the use of ERA5-Land reanalysis as precipitation input to the model versus local high-quality measurements at MF. Overall, scores at MF are deemed acceptable given that snow height is highly heterogeneous spatially under the canopy (Parajuli et al., 2020).

Figure 7 Comparison of observed (black) versus modeled snow height using SNOWPACK with the initial canopy module (IM; pink) in the Montmorency Forest (a) and in the Bernard River Valley (b).

As shown in Fig. 8a, SNOWPACK with IM simulates T_surf with an RMSE of 1 °C and a p_bias of 0.5 % in 2018 and from 2020 to 2023 (winters 2018–2019, 2020–2021, and 2021–2022 are presented in the Supplement). However, there is some discrepancy between the model and the observations when looking at a daily timescale (Fig. 8b). In this example, SNOWPACK accurately estimates T_surf from 7 to 10 March, corresponding to cloudy conditions. However, in the following 3 d, corresponding to clear-sky conditions, the model overestimates T_surf during the day and at night. In general, the results suggest that the energy balance is adequately described by SNOWPACK, so we can trust the timing of the onset of melt. For example, the observed rise to 0 °C in early April is well reproduced by the model (Fig. 8a). Figure S4 shows that this is also the case for the other years.

Figure 8 Comparison of observed (black) versus modeled snow surface temperature using SNOWPACK with the initial canopy module (IM; pink) during the winter of 2022–2023 (a) at Montmorency Forest (MF). (b) Closer look at the observed and simulated surface temperatures during the week of 7 to 14 March 2023.

The model shows an average agreement score of 0.79 for snow density profiles compared to observations at MF (Fig. 9a). Density agreement scores are generally constant throughout the season, except during winter 2018–2019 at MF, where the performance decreases during the melt period when the density is more uniform vertically. The model performs better at simulating the upper snow layers than the lower layers, where SNOWPACK tends to overestimate the density (Figs. S5 to S11). This is not surprising given that SNOWPACK does not account for upward convective water vapour fluxes (Domine et al., 2019) or vapour diffusion in the version that we used (Jafari et al., 2020). At BRV, the agreement score for snow density is 0.77 and the vertical density profile is also better represented at the top than at the bottom (Figs. S12 to S15). Although we could evaluate the model's performance based on only on six profile observations at BRV (Fig. 9b), these results are similar to those from MF and support the conclusion reached for this site.

SNOWPACK does not simulate grain type as well as snow density, with an average agreement score of 0.59 at MF (Fig. 9c). Lower agreement scores for grain type are also obtained for simulations at BRV (Fig. 9d). As with density measurements, there is no clear difference in performance between seasons for the grain type. In winter 2018–2019, there is a strong decrease in agreement score from the beginning of January to the end of April. This decrease in performance occurs immediately after an ROS event on 21 December 2018. SNOWPACK simulates the complete wetting of the snow column, whereas a thinner melt–freeze layer was observed (MF18–19b in Table 3). The score then increases during the 2019 snowmelt, as the melt forms observed in the field were also mostly simulated by the model. This is not a surprise given that melt forms are more uniform than other grain types (Colbeck, 1982).

Figure 9 Agreement scores between the SNOWPACK simulations using the initial canopy module (IM) and the observations for the snow density profiles (a, b) and snowpack stratigraphy (c, d). Agreement scores are computed for all snow pit measurements taken between 2018 and 2023 at Montmorency Forest (MF) and between 2019 and 2023 at Bernard River Valley (BRV).

4.3 Age-based intercepted snow densification

4.3.1 Effects on snow unloading

The correspondence between snow in trees inferred from photographs and the variations in modeled snow canopy storage (Fig. 10a) suggests that the removal of intercepted snow by unloading and evaporation is well captured by the model. The model reproduces well the observations of the canopy becoming snow-free due to the large drop in interception capacity triggered by ROS events, leading to large unloading events. However, ISD causes the unloaded snow to have a significantly larger density compared to IM (Fig. 10b). Note that the decrease in intercepted snow density is due to the interception of new snow, which has generally a lower density than the snow already stored in the canopy.

Figure 10 Canopy storage and intercepted snow density as simulated by SNOWPACK at Montmorency Forest (MF) in 2018–2019 with the initial canopy module (IM; pink) and the intercepted snow densification (ISD; blue; a–b). Triangles in (b) show the simulated snow density as it unloads from the canopy. The gray shaded areas in (a) and (b) indicate the presence of snow in the canopy as observed by the time lapse cameras. No photos were taken at the beginning of the season, indicated by the yellow shaded area. (c) Snow density profiles simulated with both versions of SNOWPACK and compared to observations on 5 February 2019. The simulated profiles are aggregated to a 4 cm vertical resolution and stretched to match the observed snow height.

4.3.2 Effects on snowpack properties

As shown in the density profile comparison (Fig. 10c), ISD produces near-surface density values in better agreement with observations than those of IM, which strongly underestimates density resulting from snow unloading. A higher initial density with ISD is also consequently compensated by lower settling rates. Overall, except for layers resulting from snow unloading, IM and ISD simulate similar density profiles (Figs. S5 to S15).

The use of SNOWPACK with ISD has a greater impact on the simulated grain type as shown in Fig. 11. From early December 2018 to mid-April 2019, we observe, from bottom to top, apart from melt–freeze layers, depth hoar, faceted crystals, rounded grains, and recent snow. Most importantly, we observe four melt–freeze layers resulting from ROS events (MF18–19a–d in Table 3). With IM, only the basal melt–freeze layer (MF18–19a) is well simulated, with the upper melt–freeze layers either too thick (MF18–19b), absent (MF18–19c), or not generated at the correct height (MF18–19d). In contrast, ISD reproduced all observed melt–freeze layers quite well until late March (Fig. 11c). Other thin melt–freeze layers were formed in April, but these were difficult to identify from snow pit observations and were therefore not used to evaluate model performance. Of the 18 melt–freeze layers observed at MF (Table 3), ISD is able to reproduce 17 of them, compared to 10 by IM (Figs. S16 to S19). At BRV, ISD simulated 10 layers instead of 8 by IM, supporting our results from MF (Figs. S20 to S23). All ROS-induced melt–freeze layers correctly reproduced with ISD from October through March were preceded by snow unloading. However, not all simulated unloading events were followed by an ROS event. We did not note instances where these buried layers impeded water percolation from ROS later in the season as snow metamorphism reduces the fine-over-coarse layer transitions.

Since the melt–freeze layers are a small fraction of the snow column, this improvement does not translate into a large improvement in the grain type agreement scores (Table 4), which also indicates a limitation inherent to the formulation of the agreement score. For snow height, T_surf, and density profiles, both versions of the canopy module also give similar results (Table 4). This shows that the performance gain of ISD for simulating snow stratigraphy is achieved without compromising the performance with respect to other snow properties.

Figure 11 Sub-canopy snowpack stratigraphy in the Montmorency Forest (MF) during winter 2018–2019 as (a) observed from 26 snow pits and simulated with SNOWPACK using (b) the initial canopy module (IM) and (c) the intercepted snow densification module (ISD). The following grain types are present: PP (precipitation particles), DF (decomposed and fragmented PP), RG (rounded grains), FC (faceted crystals), DH (depth hoar), SH (surface hoar), MF (melt forms), MFcr (melt–freeze layers), and IL (ice layers). MFcr and IL are shown in darker red and cyan, respectively. The colour code for snow grain type is taken from the International Classification for Seasonal Snow on the Ground (Fierz et al., 2009).

Table 4 Summary of the performance of SNOWPACK for simulating the sub-canopy snowpack in the Montmorency Forest (MF) and in the Bernard River Valley (BRV) for both versions of the canopy module (IM and ISD). The summary includes the evaluation of snow height and snow surface temperature (T_surf; RMSE, p_bias), density and grain type profiles (agreement score), and the number of ROS-induced melt–freeze layers detected. All metrics cover the entire study at both sites, except T_surf which excludes winter 2019–2020.

^a RMSE in cm for snow height and in °C for the snow surface temperature; ^b out of 18 layers at MF and 14 at BRV.

4.3.3 Case of the ROS events on 5 and 8 February 2019

To assess the timing of the effects of intercepted snow densification in simulations, we carefully examine two ROS events that took place on 5 and 8 February 2019 at the MF site (see Fig. 12). The snow that unloaded at the beginning of both ROS events produced a thick, low-density snow layer when simulated with IM (Fig. 12a). This resulted in liquid water that quickly percolated through the snowpack by preferential flow (Fig. 12b). Since water remained in the preferential flow domain, the snow column temperature stayed below 0 °C and no melt–freeze layer was created (Fig. 12c–d). As preferential flow conveys water downward rapidly, the model generated runoff a few hours only after the onset of rainfall for both ROS events (Fig. 12e). However, runoff volumes were smaller than rainwater volumes, indicating that some of precipitation was retained in the snowpack.

With the ISD module, unloaded snow triggered by rainfall created a dense layer of snow on top of the snowpack (Fig. 12f), creating fine-over-coarse conditions that restricted downward water flow. This caused water to accumulate in the preferential flow domain (Fig. 12g) until a small fraction of it percolated through when the pressure head exceeded the water entry pressure of the underlying layer with increasing grain size. The other fraction of water was transferred to the matrix flow domain when θ_TH was exceeded (Fig. 12h), which led to the formation of a melt–freeze layer (Fig. 12i). This is in line with the observed stratigraphy from 5 February 2019 presented on the right of Fig. 12i. Since water was mostly retained in the snowpack to further freeze, the ROS events resulted in almost no runoff (Fig. 12j).

Figure 12 The 5 February 2019 rain-on-snow (ROS) event in the Montmorency Forest (MF) simulated with SNOWPACK using the initial canopy module (IM; left column) and the intercepted snow densification module (ISD; right column). Panels (a) and (f) show liquid precipitation (P) on top and the evolution of snow density profile during the ROS event. Panels (b) and (g) show the evolution of volumetric liquid water content (VWC) in the preferential flow domain. Panels (c) and (h) present the evolution of snow temperature and the presence of VWC in the matrix flow domain in shaded white (also highlighted by black rectangles). Panels (d) and (i) show the snowpack stratigraphy during the ROS event, along with grain type observations from 5 February 2019. The black rectangles mark the simulated profile on this date. Finally, panels (e) and (j) show the snowpack runoff resulting from the ROS events.

4.3.4 Effects on snowpack runoff

The version of the canopy module used (IM versus ISD) dictates the simulated runoff resulting from ROS events from November through March (Fig. 13a). Compared to the simulations with IM, 70 % of the events simulated with ISD were smaller and generated 27 % less runoff volume. Of the 47 events, only 12 with IM and 9 with ISD have runoff greater than rainfall accumulation. This means that most ROS events during the accumulation period result in rainwater freezing in the snowpack rather than ROS events triggering snowmelt. For two major rainfall events at MF, the runoff simulated by ISD was 54 mm (21 December 2018) and 44 mm (24 December 2020), less than that simulated by IM. In both cases, IM generated snowmelt, while ISD caused rainwater to freeze in the snowpack, demonstrating that the canopy module has an impact on the simulated runoff. Most importantly, this suggests that canopy snow metamorphism may influence runoff generation during rain-on-snow events. Simulations in forest gaps and observations of snowpack runoff would be needed to confirm this hypothesis.

Figure 13 Simulated runoff from the initial canopy module (IM) and the intercepted snow densification module (ISD) for all 47 rain-on-snow events from November to March at both sites. The black line shows the 1:1 relationship.

4.4 Sensitivity analysis

Of the 26 snow pits of winter 2018–2019 at MF, melt–freeze layers were reported 81 times by the observer (Table 5). The number of simulated melt–freeze layers is not very sensitive to canopy snow parameters with hits ranging from 57 (ρ_max=200 kg m⁻³) to 65 (ρ_th=100 kg m⁻³). However, the number of melt–freeze layers simulated with the correct thickness is more sensitive, with hits ranging from 21 (d_g=0.1 mm) to 39 (ρ_th=200 kg m⁻³). Compared to IM, ISD detects more layers and simulates the thickness better, regardless of the sensitivity analysis parameters that we used. This suggests that ISD is more capable of simulating melt–freeze layers with the correct thickness than IM.

As shown in Fig. 14, simulated runoff appears to be more sensitive to a change in canopy snow parameters than the number of modeled melt–freeze layers. Except for one point corresponding to an event in late March, both ρ_fr and ρ_th have a very small effect on the average cumulative runoff from an ROS event (±0.3 mm). In contrast, simulations performed with a reduced and increased ρ_max result in event-based cumulative runoff that is on average 2.5 mm lower and 2.1 mm higher, respectively, than for simulations with all parameters unchanged. Similarly, decreasing and increasing d_g result in simulated cumulative runoff that is 2.9 mm higher and 4.6 mm lower, respectively, on average. This indicates that the average grain size when snow unloads has a greater effect on ROS-induced runoff than the intercepted snow density.

Table 5 Number of observed melt–freeze layers reproduced (“hits”) and simulated with the right thickness by SNOWPACK with IM, ISD, and ISD with the modified canopy parameters identified (see Table 2) in winter 2018–2019 at MF. Note that 81 melt–freeze layers were observed in snow pits.

Figure 14 Comparison of ROS-induced runoff simulated by SNOWPACK with ISD when the canopy snow parameters are unchanged (x axis) and modified (y axis). The modified parameter is identified at the top of each plot. The black line in all plots shows the 1:1 relationship.

5 Discussion

5.1 Influence of ROS events on the sub-canopy snowpack

The wet and cold conditions that prevail in the boreal forest of eastern Canada cause the snowpack to form early in the season, thicken, and persist until late May or early June. Most ROS events occurred in November–December, as temperatures from January to March were generally too low for ROS events to take place. This contrasts with observations from the western United States and central Europe, where milder conditions are conducive to frequent ROS events from January to March for areas at elevations similar to our sites (McCabe et al., 2007; Hotovy et al., 2023).

The thickness of the melt–freeze layers from ROS events shows a high variability (Table 3). This could be attributed to the amount of precipitation and air temperature during an ROS event, as well as the unloading of snow clumps that are further redistributed unevenly across the snow surface. In addition, the thickness of melt–freeze layers can decrease or even disappear over time due to gradient metamorphism (Domine et al., 2009) or increase as rainwater from subsequent ROS events accretes on the layer (Kapil et al., 2010). Overall, this complements the findings of previous studies, i.e., that vegetation modulates the factors driving ROS runoff such as the snow water equivalent and the available energy for melting the snow (Storck et al., 2002; Wayand et al., 2015).

The cases presented in Fig. 6 suggest that preferential flow is an important transport mechanism for liquid water in the sub-canopy snowpack. These results are supported by the photos in Fig. 5, which demonstrate that ROS events may result in preferential flow under the canopy. The presence of percolation channels in the sub-canopy snowpack was also observed in other studies (Bründl et al., 1999; Teich et al., 2019). Percolation channels are formed in cold snowpacks before matrix flow takes over as the snowpack warms and wet snow metamorphism prevails (Hirashima et al., 2019). Although percolation channels can form at any place in the snowpack due to boundary conditions (Schneebeli, 1995; Avanzi et al., 2016), it is more likely to find larger channels lower in the snowpack, as in Fig. 5c. This is attributed to larger grains such as faceted crystals and depth hoar that reduce the water-entry capillary pressure head of wetted snow layers (Waldner et al., 2004; Katsushima et al., 2013). In general, water percolating through the snowpack by preferential flow accelerates the hydrological response to ROS events when compared to matrix flow (Singh et al., 1997; Waldner et al., 2004; Würzer et al., 2016). Since percolation channels were widely observed under the canopy, it suggests that snow conditions and properties under the canopy favour this water transport mechanism during ROS events.

5.2 SNOWPACK simulations under a boreal canopy

Our results show that SNOWPACK generally overestimates the snow height below the canopy at both sites (Fig. 7). At BRV, this could be attributed to the accuracy of ERA5-Land precipitation used as an input variable. The overestimation of snow height could also be explained by too little interception simulated by SNOWPACK due to an underestimation of evaporation and sublimation by the model. This causes the storage capacity of the canopy to be reached too rapidly, thereby increasing snow accumulation on the ground. From 2018 to 2023 at MF, we simulated evaporation and sublimation losses ranging from 32 and 55 mm during the snow cover period. This corresponds to roughly 30 % of what Isabelle et al. (2020) observed at the same site for the snow seasons of 2016–2017 and 2017–2018. We thus conclude that SNOWPACK underestimates mass loss of intercepted snow by as much as a factor of 3, highlighting the need for further model development in this area.

The surface temperature simulated by SNOWPACK is overestimated during the day and underestimated during the night under clear-sky conditions. This could be explained by too much and too little downwelling longwave radiation under the canopy during the day and night, respectively, contributing to snow surface warming, as noted by Gouttevin et al. (2015) for simulations at a subalpine site in Switzerland with an LAI of 3.9 m² m⁻². Under cloudy conditions, SNOWPACK simulations of T_surf improve because air and snow surface temperatures are similar. Our scores for snow density and grain types are inferior to those of simulations performed in alpine terrain (Lehning et al., 2001) or forest openings (Rasmus et al., 2007). This is partly explained by difficulties in sampling snow density and identifying grain types under trees due to the heterogeneous layering of the sub-canopy snowpack (Bouchard et al., 2022). Also, two-dimensional snow–forest processes that affect snowpack structure such as snow unloading and preferential canopy dripping, are difficult to capture in a one-dimensional model.

5.3 Canopy parameterization, snowpack structure, and runoff

Using our intercepted snow densification scheme instead of the original module is a net gain in simulating the number and thickness of melt–freeze layers in the sub-canopy snowpack. Indeed, the unloading of denser snow consisting of small rounded grains creates fine-over-coarse transitions where liquid water is retained (Wever et al., 2016). In nature, melt–freeze layers are more likely promoted by a strong sub-canopy snowpack heterogeneity resulting from non-uniform processes such as unloading and meltwater dripping (Teich et al., 2019; Bouchard et al., 2022). However, our study suggests that simulating the canopy snow evolution helps reproduce the general features of melt–freeze formation under canopy (Table 5).

Delayed and reduced ROS runoff is an indirect consequence of simulating the unloading of denser snow of small rounded grains with ISD. The greater effect of grain size on runoff sensitivity than snow density (Fig. 14) can be explained by the parameterization of the water entry pressure of the snow layers, which drives the transition between the preferential and matrix flow domains (Wever et al., 2016). Since the hydraulic conductivity and the Van Genuchten parameters are estimated from ρ_s and d_g in SNOWPACK (Wever et al., 2014, 2015), snowpack runoff is also affected by the density of unloaded snow. This sheds light on the hydrological influence of microstructural descriptors of intercepted and unloaded snow. In the absence of snowpack runoff measurements, we cannot directly validate the ISD parameterization and whether this constitutes an improvement in the simulated snowpack hydrological response to ROS events. The MF site is located within a catchment in which the daily average discharge is continuously gauged at the outlet (station 051004 Des Aulnaies, https://www.cehq.gouv.qc.ca/hydrometrie/historique_donnees/default.asp, last access: 19 November 2023). Although the discharge station recorded an increase in streamflow from large ROS events of the study period, we refrain from drawing conclusions from these measurements for two reasons. First, the daily time step of discharge measurements is too coarse a resolution, as ROS-induced runoff was usually generated within 1 d of the rainfall event. Second, the forest cover in the catchment is discontinuous and snow properties are highly heterogeneous spatially, so water transport mechanisms in the snowpack may differ significantly under canopy and in forest gaps, as shown by Bouchard et al. (2022). Also, the relationship between snowpack runoff and discharge in streams should consider initial soil moisture and the flow network (Wever et al., 2017).

Nevertheless, we compared our simulations with sprinkling water experiments in subalpine snowpacks by Singh et al. (1997) and Juras et al. (2017), who showed a runoff response to rainfall within an hour. This is faster than what we observe from temperature profile measurements (Fig. 6) or simulate with ISD (Fig. 12j) with a runoff response of several hours. Conditions in those experiments such as a high rainfall rate, a shallow and warm snowpack, and the absence of canopy cover could explain the differences with our study.

Our modeling results are specific to the unloading parameterization used in SNOWPACK (see Eq. A5), which promotes large snow unloading during ROS events. Lumbrazo et al. (2022) showed that the choice of unloading parameterization affects the timing of snow unloading, as well as the canopy mass removal from unloading and sublimation. It is likely that the unloading parameterization would also influence the snow microstructure and ROS-induced runoff simulations.

5.4 Limitations

The first limitation of this work is the lack of measurements of intercepted snow properties in the canopy. Such measurements would have helped to accurately parameterize the age-based densification function. The sensitivity analysis shows relatively low impact of canopy snow parameters on melt–freeze layer formation but a larger influence on runoff that cannot be neglected (Fig. 14). In the past, observational studies have mostly focused on quantifying the mass of snow intercepted by the canopy (Friesen et al., 2015; Raleigh et al., 2022). Given the effect of canopy snow properties on ROS runoff, future studies should attempt to measured fundamental snow properties in the trees like snow density and SSA and monitor canopy snow temperature to better understand snow metamorphism in the canopy.

A second important shortcoming is the lack of continuous snowpack runoff measurements, which limits our conclusions regarding the hydrological impact of ROS. Although soil–snow interface temperature provides information on the presence of liquid water at the base of the snowpack, it does not quantify the volume of water outflow. Lysimetric measurements should be used to that end despite technical and logistical challenges (Kattelmann, 2000; Floyd and Weiler, 2008; Webb et al., 2018b). Such measurements would allow us to evaluate the performance of canopy snow densification on snowpack runoff.

This study is also limited by the single-point simulations of SNOWPACK, which contrast with the highly spatial heterogeneous character of forest snow (Bouchard et al., 2022). SNOWPACK, as any other multi-layer snow models to date, is not yet designed to reproduce the spatial heterogeneity of the snow cover from tree to tree. Therefore, further modeling developments are needed to better represent spatially variable vegetation–snow processes like interception and unloading (Vincent et al., 2018) or even radiation transfer (Jonas et al., 2020) in multi-layer, microstructure-resolving snow models. This can be achieved by coupling these detailed snow models with models that resolve tree-scale processes. Recent work by Mazzotti et al. (2023), who coupled FSM2 (Mazzotti et al., 2020) to CROCUS (Vionnet et al., 2012), looks promising to this end.

6 Conclusion

In this work, we first aimed to better document how rain-on-snow (ROS) events influence the sub-canopy snowpack structure through observations. To do so, we recorded nearly 50 ROS events, monitored the snow thermal regime, and assessed snow properties and structure from snow pit measurements over multiple winters at two boreal sites representative of different bioclimatic areas of eastern Canada. We show that rain falling through the canopy before reaching the snow cover leads to the formation of thick melt–freeze layers and that preferential flow is a major water transport mechanism in the sub-canopy snowpack.

Our second objective was to evaluate the multi-layer one-dimensional snow model SNOWPACK under a boreal canopy. Although designed to simulate alpine snow covers, SNOWPACK was found to be suitable for snowy and cold boreal environments such as those found in eastern Canada. It provides acceptable simulations of snow height and snow density and accurately simulates the snow surface temperature. However, the observed melt–freeze layers resulting from ROS events were generally not well simulated by the model.

The third objective of this study was therefore to assess the impact of implementing an age-based intercepted snow densification function in SNOWPACK on snowpack structure and runoff from ROS events. We obtained improved simulations of the number and thickness of melt–freeze layers. These improvements were obtained at both sites, illustrating the transferability of the function. The densification function also reduces and delays ROS-induced snowpack runoff as it produces dense and fine-grain layers of unloaded snow over lighter snow layers of large grains. The parameters of this function have a low impact on the simulation of melt–freeze layers. In contrast, these parameters were found to affect more strongly the simulated ROS-induced runoff, highlighting the need for documenting the physical properties of snow in the canopy.

In summary, our findings show that the evergreen canopy modulates snowpack structure, preferential flow, and snowpack runoff during ROS events. Our work is another step towards better reproducing canopy snow properties of and provides insights for further observational and modeling efforts in hydrology applied to snow-dominated forested environments. Investigating the effect of canopy snow properties on runoff at larger scales, developing a robust methodology to assess canopy snow metamorphism from observations, and coupling detailed canopy structure schemes to multi-layer snow models to simulate snow microstructure in forest gaps would be logical next steps. The approach recently developed by Mazzotti et al. (2023) looks promising for this purpose. Finally, the multi-year dataset presented in this study can further be used for future model validation and improvement in a context of increasing winter rainfall events.

Appendix A: Measurements height of the forcing variables

Table A1 Measurement height from the ground of meteorological variables at Montmorency Forest (MF) and Bernard River Valley (BRV).

Appendix B: SNOWPACK initialization file parameters

Table B1 SNOWPACK.ini file parameters.

Appendix C: SNOWPACK parameterization of interception and unloading

The temporal change in canopy storage ($\mathrm{d}I/\mathrm{d}t$; mm d⁻¹) is a function of interception rate (ΔI; mm d⁻¹), evaporation and sublimation of intercepted snow (E_int mm d⁻¹), and liquid or solid unloading from the canopy (U; mm d⁻¹):

$$\begin{array}{}\text{(C1)}& {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}}=\mathrm{\Delta }I-E_{\mathrm{int}}-U.\end{array}$$

The interception rate is calculated following Hedstrom and Pomeroy (1998):

$$\begin{array}{}\text{(C2)}& \mathrm{\Delta }I=c_{\mathrm{u}}(I_{\max }-I)\left(\mathrm{1}-e^{-\left(\frac{\left(\mathrm{1}-c_{\mathrm{f}}\right)P}{I_{\max }}\right)}\right),\end{array}$$

where c_u (–) is an unloading coefficient, set to 0.7 as suggested by Pomeroy et al. (1998); I is the mass initially stored in the canopy (mm); c_f (–) is the direct throughfall fraction which takes a value between 0 and 1; and P (mm d⁻¹) is the precipitation. I_max (mm) is the maximum capacity of the canopy (Schmidt and Gluns, 1991). When canopy interception is in the liquid phase, I_max takes the constant value of 0.3 LAI mm m⁻² as suggested by Gouttevin et al. (2015). Otherwise, I_max is estimated from the LAI (m² m⁻²), the density of new snow (ρ_new; kg m⁻³), and a tree-species-dependent interception capacity factor, i_max (mm):

$$\begin{array}{}\text{(C3)}& I_{\max }=i_{\max }\left(\mathrm{0.27}+{\displaystyle \frac{\mathrm{46}}{{\mathit{\rho }}_{\mathrm{new}}}}\right)\text{LAI},\end{array}$$

where ρ_new is a polynomial function of T_a in °C, RH (0–1), WS in m s⁻¹, and T_surf in °C:

$$\begin{array}{}\text{(C4)}& \begin{array}{rl}{\mathit{\rho }}_{\mathrm{new}}& =\mathit{\alpha }+\mathit{\beta }{T}_{\mathrm{a}}+\mathit{\gamma }{T}_{\mathrm{surf}}+\mathit{\delta }\text{RH}+\mathit{\eta }\text{WS}+\mathit{\phi }{T}_{\mathrm{a}}{T}_{\mathrm{surf}}\\ & +\mathit{\mu }{T}_{\mathrm{a}}\text{WS}+\mathit{\nu }\text{RHWS}+o{T}_{\mathrm{a}}{T}_{\mathrm{surf}}\text{RH},\end{array}\end{array}$$

where α=90, β=6.5, γ=7.5, δ=0.26, η=13, $\mathit{\phi }=-\mathrm{4.5}$, $\mathit{\mu }=-\mathrm{0.65}$, $\mathit{\nu }=-\mathrm{0.17}$, and o=0.06. ρ_s,int is limited to values between 30 and 250 kg m⁻³.

Unloading is the difference between the canopy storage and I_max over each time step and happens when the storage exceeds the maximum capacity of the canopy:

$$\begin{array}{}\text{(C5)}& U={\displaystyle \frac{max\left[\mathrm{0},I-I_{\max }\right]}{\mathrm{\Delta }t}}.\end{array}$$

The throughfall is calculated as follows:

$$\begin{array}{}\text{(C6)}& T=P-\mathrm{\Delta }I+U.\end{array}$$

In this parameterization, unloading and precipitation are merged before being added as a new snow layer to the snowpack. Unloading from the canopy takes therefore the same properties as solid or liquid precipitation.