Multilayer observation and estimation of the snowpack cold content in a humid boreal coniferous forest in eastern Canada

Cold content (CC) is an internal energy state within a snowpack and is defined by the energy deficit required to attain isothermal snowmelt temperature (0°C). For any snowpack, fulfilling the cold content deficit is a pre-requisite before the onset of the snowmelt. Cold content for a given snowpack thus plays a critical role because it affects both the timing and the rate of snowmelt. Estimating the cold content is a labour-intensive task as it requires extracting in-situ snow temperature 10 and density. Hence, few studies have focused on characterizing this snowpack variable. This study describes the multilayer cold content of a snowpack and its variability across four sites with contrasting canopy structures within a coniferous boreal forest in southern Québec, Canada, throughout winter 2017-18. The analysis was divided into two steps. In the first step, the observed CC data from weekly snowpits for 60% of the snow cover period were examined. During the second step, a reconstructed time series of CC was produced and analyzed to highlight the high-resolution temporal variability of CC for the 15 full snow cover period. To accomplish this, the Canadian Land Surface Scheme (CLASS; featuring a single-layer snow model) was first implemented to obtain simulations of the average snow density at each of the four sites. Next, an empirical procedure was used to produce realistic density profiles, which, when combined with in situ continuous snow temperature measurements from an automatic profiling station, provides a time series of CC estimates at half-hour intervals for the entire winter. At the four sites, snow persisted on the ground for 218 days, with melt events occurring on 42 of those days. Based on snowpit 20 observations, the largest mean CC (−2.62 MJ m) was observed at the site with the thickest snow cover. The maximum difference in mean CC between the four study sites was −0.47 MJ m, representing a site-to-site variability of 20%. Before analyzing the reconstructed CC time series, a comparison with snowpit data confirmed that CLASS yielded reasonable estimates of the snow water equivalent (SWE) (R = 0.64 and percent bias (Pbias) = −17.1%), bulk snow density (R = 0.71 and Pbias =1.6%), and bulk cold content (R = 0.90 and Pbias = −2.0%). A snow density profile derived by utilizing an 25 empirical formulation also provided reasonable estimates of cold content (R = 0.42 and Pbias = 5.17%). Thanks to these encouraging results, the reconstructed and continuous CC series could be analyzed at the four sites, revealing the impact of rain-on-snow and cold air pooling episodes on the variation of CC. The continuous multilayer cold content time series also provided us with information about the effect of stand structure, local topography, and meteorological conditions on cold content variability. Additionally, a weak relationship between canopy structure and CC was identified. 30


Introduction
The use of spatially distributed, process-based (physical) hydrological models has substantially improved decision-making in the area of water resources management (Wigmosta et al., 2002). The snow processes included in such models rely on the energy balance (EB) approach, since snow accumulation and melt depend on the exchanges of energy and mass between the 35 snowpack and its surrounding environment (soil, atmosphere, and vegetation). The concept of snowpack energy budget was first introduced by the U.S. Army Corps of Engineers (1956). Since then, the single bulk layer representation (e.g Wigmosta et al., 1994) has evolved into multilayer schemes (Gouttevin et al., 2015;Koivusalo et al., 2001;Lehning et al., 2002;Vionnet et al., 2012). Recent studies have looked at the sources of uncertainty associated with snow models (Essery et al., 2013;Rutter et al., 2009) and revealed the importance of including some key state variables, particularly cold content, in their modelling 40 schemes.
Cold content (CC) is the amount of energy required for a snow cover to reach 0°C for its entire depth. Any additional energy input translates into melting. By definition, CC is a linear function of the snow water equivalent (SWE) and snowpack temperature, and is defined by: where CC is cold content (MJ m −2 ), ci is the specific heat of ice (2.1 × 10 −3 MJ kg −1 °C −1 ), ρs is the snow density (kg m −3 ), ρw is the density of water (kg m −3 ), HS is the snow depth (m), Tm is the melting temperature (0°C), and Ts is the snowpack temperature (°C). Thus, CC ranges from -∞ to 0 MJ m −2 , meaning that the larger the absolute value of CC, the more energy required for the snowpack to eventually reach a uniform temperature of 0°C, over the entire depth of the snowpack.
CC plays a central role in delaying snowmelt (Molotch et al., 2009), as a deep, dense, and cold snowpack requires a substantial 50 amount of energy for snow to reach 0°C and initiate melt. As such, understanding CC is essential for the accurate forecasting of water availability in demanding sectors such as agricultural systems, urban water supply (Barnett et al., 2005), and hydropower generation (Schaefli et al., 2007).
The exact determination of CC requires direct observations of the snowpack temperature, density, and depth, usually collected from manual snow surveys. As manual collection is tedious and demanding, few datasets that describe snowpack CC are 55 available. For lack of a better approach, CC is often estimated using one of the following three methods: an empirical formulation that relies solely on air temperature (DeWalle and Rango, 2008;Seligman et al., 2014), an empirical formulation based on air temperature and precipitation (Andreadis et al., 2009;Wigmosta et al., 1994), or a residual from an energy balance model (Marks and Winstral, 2001). Jennings et al. (2018) resorted to snowpit data, collected at alpine and subalpine sites within the Rocky Mountains in Colorado, to study CC. They reported a weak relationship between CC and the cumulative 60 mean of air temperature. They also found that newly fallen snow was responsible for 84.4% and 73.0% of the daily gains in CC for alpine and subalpine snowpacks, respectively. Of note, a slight contrast was observed by Seligman et al. (2014), who reported that the contribution of spring snow storms to CC had a smaller impact on delaying snowmelt than the porous space from dry fresh snow. However, Jennings et al. (2018) reported shifts in the onset of snowmelt by 5.7 h and 6.7 h at alpine and https://doi.org/10.5194/tc-2021-98 Preprint. Discussion started: 13 April 2021 c Author(s) 2021. CC BY 4.0 License. subalpine sites, respectively, when CC at 6:00AM was less than 0 MJ m −2 . This suggests that even a small energy deficit has 65 a substantial effect on the rate and timing of snowmelt. Overall, previous studies agree that the careful consideration of CC improves snowmelt simulations (Jost et al., 2012;Mosier et al., 2016;Valéry et al., 2014).
Little to no previous research has focused on CC behaviour in forested environments. Snowpack energy exchanges within a forest are obviously different than those in open or alpine areas, as the presence of a canopy impacts snow accumulation and melt (Andreadis et al., 2009;Gouttevin et al., 2015;Mahat and Tarboton, 2012;Wigmosta et al., 2002). For instance, 70 intercepted snow may sublimate, undergo densification, or fall beneath the canopy when maximum canopy storage is reached or when there are heavy winds present (DeWalle and Rango, 2008). Snow interception typically leads to shallower snow depths and less melt beneath the canopy (Musselman et al., 2008), even in the presence of rain-on-snow events (Marks et al., 1998).
Frequent density profiles of the snow cover allow for the tracking of unloading episodes and the identification of spatial differences of CC within a forest. 75 Despite all of the associated challenges, it is possible to simulate snow in a forested environment with some success. For instance, physically-based land surface models are regularly used to simulate snow at forested sites (e.g., Roy et al., 2013).
One such example is the Canadian LAnd Surface Scheme (CLASS), which relies on a single-layer snow model articulated around the energy balance. In a recent study, Alves et al.(2020) used CLASS driven by ERA5 reanalysis data to model snow depths from four dissimilar forested sites across the Canadian boreal biome. They reported average snow persistence lengths 80 and average spring melting periods that were similar to our field observations. By definition, CLASS considers the whole snowpack as a single bulk unit, and as such, is unable to simulate the multilayer behaviours that one sees in nature. One option for addressing this is to resort to a multilayer snow model such as SNOWPACK (Lehning et al. 2001), which was recently equipped with a thorough canopy module (Gouttevin et al., 2015); however, even models such as this are not free of biases.
Alternatively, bulk snowpack values can be distributed between several layers. For instance, Roy et al. (2013) disaggregated 85 CLASS-derived snow water equivalents into multilayer values at each time step, for the purpose of estimating the specific surface area (SSA) of a snowpack. They attained an acceptable root mean square error (RMSE) of 8.0 m 2 kg −1 in CLASSderived SSA for individual layers.
In view of the obvious lack of observational studies that are required to support model development in forested environments, detailed analyses of multilayer in situ snowpack CC are necessary. Building on Jennings et al. (2018), this study investigates 90 53 snowpit-derived CC observations at four distinct coniferous forested sites, over the course of one winter. The temporal variability of the CC is also analysed by reconstructing time series that include bulk and multilayer CC with a 30-min time step, and combine automated snow temperature observations and bulk snow density estimates that was calculated using the CLASS model.

Study sites and data collection
Observations were collected in the Bassin Expérimental du Ruisseau des Eaux-Volées (BEREV), which is a small boreal forest catchment within Montmorency Forest, Quebec, Canada (Fig. 1). This region experiences substantial precipitation (1583 mm), with 40% falling in solid form between November and May (Isabelle et al., 2018). The boreal catchment lies in the Laurentian Mountains of the Canadian Shield and is characterized by a humid continental climate (Schilling et al., 2021). There are patches 100 of forest clearings found within the basin due to past logging operations that have led to variability in stand structure (Parajuli et al., 2020b). Over the years, several vegetation species such as black spruce (Picea mariana (Mill.)) and white spruce (Picea glauca (Moench.)) were planted. However, the environment favoured the regrowth of balsam fir (Abies balsamea) stands. Isabelle et al. (2020) provide detailed information on the vegetation cover at the study site. The current analysis focuses on the four contrasting sites presented in Table 1. 110 Inspired by Lundquist and Lott (2010), we deployed an automated snow-profiling station at each location, composed of 18 Ttype thermocouples vertically spaced 10 cm apart and of an ultrasonic depth sensor (Judd Communication, USA). An additional T-type thermocouple was enclosed in a radiation shield (Fig.1c) 2 m above ground for simultaneous air temperature measurements. Maintaining a weekly timeline was sometimes difficult due to uncontrollable circumstances such as freezing rain, rain-onsnow events or even winter storms. During melt, from 21 April 2018 and on, it was impossible to reach all study sites because 120 of reduced snow depths preventing the safe use of snowmobiles, except for site A2 that was more easily accessible from the main road. Snow-profiling stations malfunctioned occasionally (less than 1% of the time), mostly in spring. Missing values were filled with snowpit observations. An exponential moving average procedure was implemented to reduce noise in snow depth observations.

Construction of CC time series 125
The exercise of constructing 30-minute time series of the snowpack CC represents a certain challenge. On the one hand, it requires time series of the vertical profile of snow temperature, which is obtained from the snow-profiling stations. On the other hand, time series of the snow density profile are needed as well. This is where the main difficulty lies. A simple approach would be to interpolate the density values extracted from snowpits, but this would be incomplete and error-prone given their limited number and absence early in the season. Herein, it was opted to produce multilayer time series of snow density thanks 130 to CLASS bulk simulations complemented with empirical formulations, as detailed below.

CLASS model
CLASS is a physically-based land surface model that simulates the exchanges of water and energy between the Earth's surface and the atmosphere (Bartlett et al., 2006;Verseghy, 1991). It considers four distinct surface subareas: bare soil, canopy cover over bare soil, canopy with snow cover, and snow cover over bare soil (Bartlett and Verseghy, 2015;Verseghy et al., 2017). 135 In this analysis, CLASS version 3.6 was used in offline mode and with a 30-min time step, enabling the stability of the prognostic modelled variables (Roy et al., 2013). CLASS allows for the inclusion of multiple soil layers and accounts for snow interception, snow thermal conductivity, and snow albedo, as described in Bartlett et al. (2006). The following subsections describe the meteorological forcing data required to run CLASS, and the methodology (CLASS + snow-profiling station) to produce single-and multi-layer time series of snow density, following Andreadis et al. (2009). 140

CLASS setup and forcing
The meteorological inputs required to run CLASS include precipitation rate, wind speed, air specific humidity, incoming shortwave and longwave radiation, air temperature, and surface atmospheric pressure (Alves et al., 2019; Leonardini et al., 2020). As CLASS is designed to explicitly consider the energy exchanges between the soil surface, vegetation, snowpack, and atmosphere, above-canopy meteorological forcings are used. The model accounts for local effects associated with the presence 145 of a canopy (e.g., attenuation of incident radiation, etc.), and incorporates user-defined parameters such as vegetation height, canopy density, and leaf area index (Table 2). Precipitation rates were determined using a GEONOR weighting gauge equipped with a single Alter shield approximately 4 km north of the study area, and were considered to be uniformly distributed throughout the catchment. Given the known windinduced bias associated with this type of gauge (Pierre et al., 2019), a simple adjustment was applied. This adjustment involved twice-daily manual precipitation observations from a Double Fence Intercomparison Reference (DFIR) setup close by, as in 155 Parajuli et al. (2020a). Vegetation parameters were extracted at each site from a LiDAR dataset. Wind speed, air specific humidity, shortwave and longwave radiation, and surface atmospheric pressure measurements were taken from flux towers at sites A1 and A2. Comparable data were unavailable at sites A3 and A4 (Table 2).
This study was carried out in a small experimental watershed with an area of 3.49 km 2 , where the sampling sites (A1 to A4) were close to one another ( Fig.1) but had distinct characteristics (Table 1). Given the similarity (more or less) in canopy 160 structure, we opted to use the inputs recorded at site A2 to run CLASS simulations at sites that lacked direct measurements of meteorological inputs (Table 2). Here we assumed negligible differences in the above-canopy inputs between sites A2, A3 and A4. The following sub-section highlights the steps adopted to generate the multilayer density estimates needed to calculate the CC time series for all snow layers.

Reconstruction of multilayer snow density time series (a hybrid approach) 165
The empirical formulation described in Andreadis et al. (2009), based on Anderson (1976), is used to reconstruct multiple layer snow density estimates by combining the CLASS-derived snow water equivalent (SWE) estimates (hereafter referred to as the hybrid procedure). Fresh snow density follows the formulation from Brun et al. (1989), who developed the method using data collected in the French Alps. We initialized the density of fresh snow by imposing a minimum snow density of 76 kg m -3 , based on available snowpit observations, and then using the equation: 170 where ρf is the density of fresh snow (kg m −3 ), um is the wind speed (m s −1 ), and Ta is the air temperature (K). With the exception of Eq 2, both Ta and Ts are presented in degrees Celsius (°C). As snow undergoes compaction due to metamorphism and the increasing weight of overlying snow, density is assumed to increase according to the following rate: where t is time (s), CRm is the snow compaction due to metamorphism (kg m -3 s -1 ), and CRo is the compaction due to the weight of overlying snow (kg m -3 s -1 ). CRm is then calculated as (Andreadis et al., 2009 where n0 = 3.6×10 −6 N s m −2 is the snow viscosity, c5 = 0.08 K −1 , c6 = 0.021 m 3 kg −1 and Ps (N s m −2 ) is the load pressure for each layer. The load pressure is defined as: where g is the acceleration due to gravity 9.8 m s −2 , Wns and Ws are the amount of new snow and the snow (derived from 185 CLASS) within the snowpack layer (mm w.e.), respectively, and f is the empirical compaction coefficient taken as 0.6 (Andreadis et al., 2009).

Local meteorological conditions
Figure 2 displays daily air temperature and wind speed observations. The shaded zone and site-specific dots illustrate the 190 temporal distribution of the manual snow surveys. Air temperature measurements were taken at 2 m above ground. Wind speed sensors were located 3 m and 2 m above ground at sites A1 and A2, respectively. To compensate for this height difference and enable fair comparisons between sites A1 and A2, wind speed measurements at site A1 were adjusted to a 2-m height, assuming a log profile. As expected, the sapling site (mean canopy of 1.8 m) experienced higher wind speeds (mean 1.3 m s -1 ) than the juvenile one (mean canopy of 8.1 m and mean wind speed of 0.12 m s -1 ). Air temperatures were homogenous from site to site, 195 with average values of −6.1 °C, −6.3 °C, −6.9 °C, and −6.5 °C at sites A1, A2, A3, and A4, respectively.   Figure 3 illustrates 10-cm CC derived from the snowpit surveys. One has to sum up all the values of a single profile to find the total CC for a specific date. Variability in snow depth, mainly induced by contrasting canopy structure, is indicated in 205 Figure 3. When comparing layer-wise (10 cm) differences (A1, A2, A3, and A4), the lowest (−0.013 MJ m −2 ) and peak (−0.67 MJ m −2 ) CC both occurred at site A2. Unlike for spring melt, when CC is low and relatively uniform, the accumulation period portrays substantial layer-wise variability structured around three distinct layers. For instance, sites A1, A2, A3, and A4 reported 9, 15, 12, and 5 observations of CC that were below −0.35 MJ m −2 . Note that the occurrence of large amplitude of CC values was not always confined to the topmost layer, as the layer just beneath the top layer also exhibited such amplitude 210 (see Fig. 3, week 10, site A3 at a snow depth of 106 cm, for example). However, the layers that are close to the ground experienced smaller amplitude of CC throughout the winter. Peak CC occurred in early February ( Fig. 4 and Table 3), when the minimum daily air temperature fell to about -25°C (Fig. 2). At that time, the amplitude of CC was highest at site A1, which also had the deepest snowpack (128 cm). The total CC time series highlighted the variability of CC across the four study sites (Fig. 4). Overall, maximum snow depth occurred at site A2 (194 cm), which also experienced the largest amplitude of mean total CC 220 (−2.62 MJ m -2 ), as a thicker snowpack can hold more CC (Table 3). For its part, A4 experienced the smallest maximum snow depth (142 cm) and the lowest amplitude of mean total CC (−2.15 MJ m −2 ). The CC difference across sites reached 0.47 MJ m -2 in total cold content, representing a variability of 20%.

Figure 5. Observed versus CLASS-simulated bulk values of SWE, snow density, and CC. R 2 denotes the coefficient of determination and Pbias (%) represents the percent bias.
After confirming that CLASS successfully models bulk snow cover variables, we moved forward with the next step in our methodolgy. We adopted a "hybrid procedure" in which we reproduced a vertical structure following Andreadis et al. (2009).

Reconstructed CC time series
Additionally, sites with less vegetation (site A1) experienced higher peak CC than sites with mature forest (A4) (Fig. 8).
Notably, the rain-on-snow episodes that occurred on 11 January, 20 February, and 30 March 2018 (thin vertical bands of low CC) were absent from the weekly series shown in Figure 3. Sites A1, A3 and A4 had a shallower snowpack and the rainfall 265 penetrated deeper, resulting in a reduced CC throughout the snowpack. Contrarily, at site A2 which had a deeper snowpack, similar rain penetration into the snowpack was only observed on 11 January 2018. All snowpacks became isothermal from 21 https://doi.org/10.5194/tc-2021-98 Preprint. Discussion started: 13 April 2021 c Author(s) 2021. CC BY 4.0 License.

Figure 8. Seasonal variability of 10-cm CC simulations stored at 30-min time intervals. The colour bar indicates CC values in MJ m −2 . Light green shading represents rain-on-snow events and light blue shading represents melt.
Due to differences in snow accumulation and melt patterns, mostly induced by differences in vegetation and topographic characteristics, there is noticeable site-to-site variability in CC (Fig. 8). The detailed variability of total CC across the four forested sites is presented in Figure 9, along with snow depth. The amplitude of total CC at site A2 was larger than at A1 275 approximately 60% of the time. At site A3, this fraction drops to 32%.

280
We examined the relationship between CC and the snow density (ρs), snow depth (HS), snowpack temperature (Ts), and air temperature (Ta; Fig. 10). Pearson's correlation coefficient (r) was used to determine these relationships for the observed and estimated values, respectively. Snowpack temperature (r = 0.83 and 0.69) and air temperature (r = 0.56, and 0.66) exhibited a positive correlation. Conversely, snow depth (r = −0.5 and −0.45) exhibited a negative correlation, whereas snow density (r = 0.4 and 0.24) showed a weak relationship. 285 Next, we examined the relationship between each of the above-mentioned variables and CC at the individual sites. This was done to identify any trends in the site-wise relationship between CC and ρs, HS, Ts, and Ta. A decreasing trend in the correlation coefficient (r) with increasing mean tree heights was observed when we examined the snow temperature and the reconstructed cold content for each site (r = 0.75, 0.69. 0.67 and 0.60 for sites A1, A2, A3, and A4, respectively). Beyond that relationship, https://doi.org/10.5194/tc-2021-98 Preprint. Discussion started: 13 April 2021 c Author(s) 2021. CC BY 4.0 License.
we did not identify any site-wise trends between CC and the other variables, thereby suggesting a weak dependency on forest 290 structure in the relationship between CC and other pertinent variables.

CC observations
As illustrated in Figures 2 and 3, the four experimental sites exhibited unique snow depths, wind speeds, and air temperatures that ultimately resulted in temporal and spatial differences in CC. Variability was such that the maximum CC was not always exhibited by the top layer, but also by the middle layer (Fig. 11). For instance, in week 15, the snowpack was denser in the top https://doi.org/10.5194/tc-2021-98 Preprint. Discussion started: 13 April 2021 c Author(s) 2021. CC BY 4.0 License. layer than in the middle layer. In week 13, the top layer snowpack was warmer than the layer beneath it. Such patterns in 300 temperature and density are counterexamples of the general patterns depicted in Figure 11. Furthermore, the importance of snow mass on CC (total) at the study sites is highlighted in Figure 4. As expected, a deeper snowpack is typically associated with higher CC. For instance, CC peaked at all sites in February, but more CC was observed 305 in the deeper A1 and A2 snowpacks. The same finding holds when CC is averaged over the 15-week period: A2 experienced more snow and higher CC, followed by A1. In both instances (peak and average CC conditions), a deeper snowpack led to larger amplitude of CC. In a similar study of alpine and subalpine snowpacks in the Rocky Mountains of Colorado, USA, Jennings et al. (2018) reported peak CC to be 2.6 times greater for the alpine snowpack than for the subalpine location, which they mostly attributed to the higher SWE accumulation at the alpine site. 310 In early February (during peak CC conditions), the snow depth difference between sites A1 and A2 was very small (Table 3).
Nonetheless, A1 exhibited higher CC than A2 (Fig. 4). This is because in addition to snow depth, CC values depend on the density and temperature of the snow (Fig. 11). The higher peak CC found at site A1 can be explained by the higher snow https://doi.org/10.5194/tc-2021-98 Preprint. Discussion started: 13 April 2021 c Author(s) 2021. CC BY 4.0 License. density that is typically associated with higher wind velocities (Vionnet et al., 2012) and wind speed-induced densification.
As illustrated in Figure 2, site A1 was windier than A2. This is expected, as it is well known that wind speed is low within 315 forest canopies (Davis et al., 1997;Harding and Pomeroy, 1996), such as those in site A2.

Reconstructed CC time-series
As mentioned previously, gaps between weekly snowpit surveys failed to capture short-lived events such as warm and cold spells or rain-on-snow events. In an attempt to produce higher frequency CC time series, we used the CLASS land surface model to simulate 30-min bulk snow density and SWE (Fig. 5). Given the limitations of bulk estimations, which are often too 320 broad to properly describe all snowpack processes (Roy et al., 2013), several studies have opted for a multilayer snow model (Brun et al., 1997;Lehning et al., 2002;Vionnet et al., 2012). Our study explored the (simpler) hybrid procedure proposed by Andreadis et al. (2009). Using this method, we generated snow density values which support the derivation of CC time series that are more prone to capturing short-lived events (Fig. 8). As reported in weekly snowpit surveys, the simulated CC time series suggest that the highest peak in CC occurred at site A1 and the highest peak in mean CC was at site A2. In both instances 325 (snowpit observations and reconstructed CC), there was a decrease in peak CC with an increase in tree height ( Fig. 8 and Table   3). Parajuli et al. (2020b) explored the spatiotemporal variability of SWE in the same forest and reported reduced snow accumulation beneath canopies composed of taller trees. Initially, site A1, with lower vegetation, experienced more snow than the other sites (Table 3). Favourable conditions (lower temperature and deeper and denser snowpack) supported the occurrence of higher peak CC at this site (Fig. 2, Fig. 11 and Table 3). However, over the entire winter, site A2 experienced more snow 330 and higher amplitude of CC values than the snowpack at site A1 (Fig. 9). In order to understand this variability in CC across all four sites, our analysis revealed that there were 60% occurrences where amplitude of CC was higher at site A2 than at A1, contributing to the overall CC at site A2.
For site A3, there were 32% occurrences where amplitude of CC values was higher than at A1, beginning in early February and continuing through the rest of the study period (Fig. 9). Most of the time, the measured snow depth at site A3 was also 335 shallower than at site A1. We hypothesized that cold air pooling might explain this phenomenon. During stable atmospheric boundary layer conditions, with weak synoptic forcing, there is reduced wind flow. This results in thermal decoupling in the valley depression, which favours the formation of a cold air pool (Fujita et al., 2010;Mott et al., 2016). This is substantiated by the rapid cooling of near-surface air within the valley depression, typically at night or early in the morning (Smith et al., 2010). As site A3 is situated in a valley depression (Fig. 1), cold air pooling most likely explains the higher peak CC at this 340 location (Fig. 2).
Based on Figures 2 and 9, snow depth and air temperature appear to influence CC distribution across the study sites. In general, the observed and simulated snowpack CC values at all sites were strongly (positively) correlated with snowpack temperature and air temperature, and weakly correlated with the snow density and snow depth values ( Figure 10). It should also be noted that the snowpack CC values at all sites only showed negative correlations with snow depth (Fig. 10). Based on CC 345 observations and the hybrid procedure, we were able to identify a relationship between the mean CC and the tree height (  Fig. 8). However, we were unable to report any trends in the site-wise relationship between CC and the above-mentioned variables (Fig. 10). Conversely, Jennings et al. (2018) attempted to establish a relationship between CC development and the cumulative mean of air temperature across the alpine and sub-alpine sites in the Rocky Mountains in Colorado, USA, but were unsuccessful. 350

Sources of uncertainty
One of the shortcomings of our multilayer snowpack scheme is the use of empirical fresh snow density estimates. Russell et al. (2020) explored a range of fresh snow density formulations and concluded that a constant value of 100 kg m −3 provided a better outcome than most empirical formulations. Nonetheless, they tested some empirical formulations that omitted the influence of wind speed on snow densification, as in the Brun et al. (1989) method. We compare observations to snow density 355 estimates (top 10 cm) from three empirical methods: Diamond-Lowry (Russell et al., 2020), Hedstrom-Pomeroy (Hedstrom and Pomeroy, 1998), and Brun (Shrestha et al., 2010;Vionnet et al., 2012) (Fig. 11). These formulations are typically used to determine snowpack density. Examples include a study by Gouttevin et al. (2015) where the Brun method was used, and one by Bartlett et al. (2006) where the Hedstrom-Pomeroy method was implemented.  (Fig. 7). Several multilayer snow models use the Brun formula/method to estimate fresh snow density (e.g. Shrestha et al., 2010;Vionnet et al., 2012). When the same (Brun method) empirical snow density model was adopted into our methodology, the snow density estimates (especially the top layer) were 365 disastrous (Figure 7). https://doi.org/10.5194/tc-2021-98 Preprint. Discussion started: 13 April 2021 c Author(s) 2021. CC BY 4.0 License.
In a slightly different context, Raleigh and Small (2017) concluded that snow density modelling was a major source of uncertainty when studying catchment SWE derived from satellite data. Additionally, several errors and biases could arise due to poor data quality and modelling deficiencies, thereby affecting the snowmelt models (Parajuli et al., 2020a;Raleigh et al., 2015Raleigh et al., , 2016Rutter et al., 2009). For instance, Jennings et al. (2018) applied the SNOWPACK multilayer model and reported 370 an overestimation in fresh-snow temperature. As reported in the present study, CC depends heavily on snowpack temperature. The quality of model inputs also influences model performance. For example, sites A1 and A2 benefitted from local flux tower measurements, but such direct measurements were not available for sites A3 and A4, for which many assumptions were necessary in order to create/complete missing input time series. This problem has also been observed in several other studies(e.g. Pomeroy et al., 2007;Qi et al., 2017). Important snowpack properties beyond just CC, such as thermal conductivity 375 (Oldroyd et al., 2013) and snow interception (Hedstrom and Pomeroy, 1998) also need to be further addressed. Therefore, future research that utilizes the physically-based snow model and describes the internal snowpack processes should focus on improving snow density estimations.

Conclusion
The purpose of this study was to document the spatial variability of CC in a humid boreal forest, using detailed measurements 380 supplemented by physically-based and empirical model outputs. The studied boreal forest is characterized by a non-uniform stand structure that led to site-to-site variations in the 10-cm weekly observations of CC. Areas with lower vegetation had the highest snow accumulation and thus resulted in the largest peaks in total CC, while the juvenile forest experienced the highest amplitude of average CC over the 15 weeks.
The Canadian Land Surface Scheme model was then coupled with complementary empirical formulations to construct bulk, 385 followed by 10-cm, 30-min snow density time series. Both CLASS and the empirical formulations supplied reasonable snow density and CC estimates. When the latter 10-cm time series were split into three layers, the bottom and the middle layers also resulted in reasonable simulations. However, modelling of the top layer was not as successful. The constructed time series were used to illustrate the influence of phenomena that are not detectable when only snowpit data are used, such as rain-onsnow episodes or the formation of cold air pools at the bottom of the valley. 390 We used the Pearson's correlation coefficient (r) to identify the role of pertinent variables (snow density, snowpack temperature, snow depth and air temperature) that affect the distribution of CC at our boreal forest sites. Snowpack and air temperature appeared to be highly influential on CC distribution compared to the depth and the density of the snowpack. Our study was supported by 30-min time step time series of 10-cm snow temperature profiles and bias-corrected precipitation inputs. The inclusion of such inputs helped us to reduce errors and biases. This study also highlighted the uncertainty associated 395 with fresh snow density estimates when simulating the physically-based snowmelt models.

Data availability:
The data that support the findings in this study will be available in the public repository.