CentrEau - Quebec Water Management Research Centre

Abstract Snowmelt contributes a significant fraction of groundwater recharge in snow‐dominated regions, making its accurate quantification crucial for sustainable water resources management. While several components of the hydrological cycle can be measured directly, catchment‐scale recharge can only be quantified indirectly. Stable water isotopes are often used as tracers to estimate snowmelt recharge, even though estimates based on stable water isotopes are biased due to the large variations of δ 2 H and δ 18 O in snow and the difficulty to measure snowmelt directly. To overcome this gap, a new tracer method based on on‐site measurements of dissolved He, 40 Ar, 84 Kr, N 2 , O 2 , and CO 2 is presented. The new method was developed alongside classical tracer methods (stable water isotopes, 222 Rn, 3 H/ 3 He) in a highly instrumented boreal catchment. By revealing (noble gas) recharge temperatures and excess air, dissolved gases allow (i) the contribution of snowmelt to recharge, (ii) the temporal recharge dynamics, and (iii) the primary recharge pathways to be identified. In contrast to stable water isotopes, which produced highly inconsistent snowmelt recharge estimates for the experimental catchment, dissolved gases produced consistent estimates even when the temperature of snowmelt during recharge was not precisely known. As dissolved gases are not controlled by the same processes as stable water isotopes, they are not prone to the same biases and represent a highly complementary tracer method for the quantification of snowmelt recharge dynamics in snow‐dominated regions. Furthermore, an observed systematic depletion of N 2 in groundwater provides new evidence for the pathways of biological N‐fixation in boreal forest soils.

Digital Twins (DTs) are on the rise as innovative, powerful technologies to harness the power of digitalisation in the WRRF sector. The lack of consensus and understanding when it comes to the definition, perceived benefits and technological needs of DTs is hampering their widespread development and application. Transitioning from traditional WRRF modelling practice into DT applications raises a number of important questions: When is a model's predictive power acceptable for a DT? Which modelling frameworks are most suited for DT applications? Which data structures are needed to efficiently feed data to a DT? How do we keep the DT up to date and relevant? Who will be the main users of DTs and how to get them involved? How do DTs push the water sector to evolve? This paper provides an overview of the state-of-the-art, challenges, good practices, development needs and transformative capacity of DTs for WRRF applications.

In snow-fed catchments, it is crucial to monitor and model the snow water equivalent (SWE), particularly when simulating the melt water runoff. SWE distribution can, however, be highly heterogeneous, particularly in forested environments. Within these locations, scant studies have explored the spatiotemporal variability in SWE in relation with vegetation characteristics, with only few successful attempts. The aim of this paper is to fill this knowledge gap, through a detailed monitoring at nine locations within a 3.49 km2 forested catchment in southern Québec, Canada (47°N, 71°W). The catchment receives an annual average of 633 mm of solid precipitation and is predominantly covered with balsam fir stands. Extracted from intensive field campaign and high-resolution LiDAR data, this study explores the effect of fine scale forest features (tree height, tree diameter, canopy density, leaf area index [LAI], tree density and gap fraction) on the spatiotemporal variability in the SWE distribution. A nested stratified random sampling design was adopted to quantify small-scale variability across the catchment and 1810 manual snow samples were collected throughout the consecutive winters of 2016–17 and 2017–18. This study explored the variability of SWE using coefficients of variation (CV) and relating to the LAI. We also present existing spatiotemporal differences in maximum snow depth across different stands and its relationship with average tree diameter. Furthermore, exploiting key vegetation characteristics, this paper explores different approaches to model SWE, such as multiple linear regression, binary regression tree and neural networks (NN). We were unable to establish any relationship between the CV of SWE and the LAI. However, we observed an increase in maximum snow depth with decreasing tree diameter, suggesting an association between these variables. NN modelling (Nash-Sutcliffe efficiency [NSE] = 0.71) revealed that, snow depth, snowpack age and forest characteristics (tree diameter and tree density) are key controlling variables on SWE. Using only variables that are deemed to be more readily available (snow depth, tree height, snowpack age and elevation), NN performance falls by only 7% (NSE = 0.66).

• Study on micropollutant distribution between biochar, APL and pyrolysis oil. • Determined PFAS in biochar, APL and pyrolysis oil. • Biochar met the guidelines for agricultural valorisation, except for heavy metals. • High methane production was achieved during anaerobic digestion of APL. • CO 2 supply during pyrolysis achieved highest calorific value of pyrolysis oil. Pyrolysis is a suitable process for sewage sludge valorisation while potentially reducing micropollutants. However, previous studies mainly focused on micropollutants in biochar, neglecting their presence of aqueous pyrolysis liquid (APL), pyrolysis oil, and gas. This study analyses the distribution of 66 micropollutants, including PFAS, dioxins and furans (PCDD/Fs), heavy metals (HMs), and polycyclic aromatic hydrocarbons (PAHs) across biochar, APL and pyrolysis oil. Additionally, the impact of a carrier gas supply (N 2 and CO 2 ) on micropollutants distribution was evaluated. In a second stage, the safe use of the pyrolysis end-products for agricultural and energetic purposes was explored. PFAS from biosolids were distributed in biochar (12–13 %), APL (6–7 %), and the oil phase (2–5 %). 63–74 % remained unaccounted for, possibly transferred to the gas phase, or decomposed during pyrolysis. PAHs were predominantly found in the pyrolysis oil, while PCDD/Fs were found in the biochar and pyrolysis oil. HMs were primarily found in biochar. PAH and PCDD/F values in the biochar met the European Biochar Certificate (EBC) guidelines. However, HMs surpassed the thresholds, suggesting either post-treatment or using biochar as a building material instead. Given that pyrolysis oil contains significant quantities of micropollutants, high-temperature combustion could serve for both micropollutant decomposition and energy reclamation. Energetic valorisation of APL assessed by biomethane potential tests, achieved a methane production yield of 432–450 NmL CH4 /g VS . Overall, a combined anaerobic digester and pyrolysis process, including a recirculation of the APL, the valorisation of pyrolysis oil for process heating, and the use of CO 2 from biogas as pyrolysis carrier gas, is suggested for further study.

Nowadays, modelling, automation and control are widely used for Water Resource Recovery Facilities (WRRF) upgrading and optimization. Influent generator (IG) models are used to provide relevant input time series for dynamic WRRF simulations used in these applications. Current IG models found in literature are calibrated on the basis of a single performance criterion, such as the mean percentage error or the root mean square error. This results in the IG being adequate on average but with a lack of representativeness of, for instance, the observed temporal variability of the dataset. However, adequately capturing influent variability may be important for certain types of WRRF optimization, e.g., reaction to peak loads, control system performance evaluation, etc. Therefore, in this study, a data-driven IG model is developed based on the long short-term memory (LSTM) recurrent neural network and is optimized by a multi-objective genetic algorithm for both mean percentage error and variability. Hence, the influent generator model is able to generate a time series with a probability distribution that better represents reality, thus giving a better influent description for WRRF design and operation. To further increase the variability of the generated time series and in this way approximate the true variability better, the model is extended with a random walk process.

Abstract. Accurate knowledge of snow depth distributions in forested regions is crucial for applications in hydrology and ecology. In such a context, understanding and assessing the effect of vegetation and topographic conditions on snow depth variability is required. In this study, the spatial distribution of snow depth in two agro-forested sites and one coniferous site in eastern Canada was analyzed for topographic and vegetation effects on snow accumulation. Spatially distributed snow depths were derived by unmanned aerial vehicle light detection and ranging (UAV lidar) surveys conducted in 2019 and 2020. Distinct patterns of snow accumulation and erosion in open areas (fields) versus adjacent forested areas were observed in lidar-derived snow depth maps at all sites. Omnidirectional semi-variogram analysis of snow depths showed the existence of a scale break distance of less than 10 m in the forested area at all three sites, whereas open areas showed comparatively larger scale break distances (i.e., 11–14 m). The effect of vegetation and topographic variables on the spatial variability in snow depths at each site was investigated with random forest models. Results show that the underlying topography and the wind redistribution of snow along forest edges govern the snow depth variability at agro-forested sites, while forest structure variability dominates snow depth variability in the coniferous environment. These results highlight the importance of including and better representing these processes in physically based models for accurate estimates of snowpack dynamics.

Abstract. 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.

Abstract This work gives an overview of the state-of-the-art in modelling of short-cut processes for nitrogen removal in mainstream wastewater treatment and presents future perspectives for directing research efforts in line with the needs of practice. The modelling status for deammonification (i.e., anammox-based) and nitrite-shunt processes is presented with its challenges and limitations. The importance of mathematical models for considering N2O emissions in the design and operation of short-cut nitrogen removal processes is considered as well. Modelling goals and potential benefits are presented and the needs for new and more advanced approaches are identified. Overall, this contribution presents how existing and future mathematical models can accelerate successful full-scale mainstream short-cut nitrogen removal applications.

Abstract. In the boreal forest of eastern Canada, winter temperatures are projected to increase substantially by 2100. This region is also expected to receive less solid precipitation, resulting in a reduction in snow cover thickness and duration. These changes are likely to affect hydrological processes such as snowmelt, the soil thermal regime, and snow metamorphism. The exact impact of future changes is difficult to pinpoint in the boreal forest, due to its complex structure and the fact that snow dynamics under the canopy are very different from those in the gaps. In this study, we assess the influence of a low-snow and warm winter on snowmelt dynamics, soil freezing, snowpack properties, and spring streamflow in a humid and discontinuous boreal catchment of eastern Canada (47.29° N, 71.17° W; ≈ 850 m a.m.s.l.) based on observations and SNOWPACK simulations. We monitored the soil and snow thermal regimes and sampled physical properties of the snowpack under the canopy and in two forest gaps during an exceptionally low-snow and warm winter, projected to occur more frequently in the future, and during a winter with conditions close to normal. We observe that snowmelt was earlier but slower, top soil layers were cooler, and gradient metamorphism was enhanced during the low-snow and warm winter. However, we observe that snowmelt duration increased in forest gaps, that soil freezing was enhanced only under the canopy, and that snow permeability increased more strongly under the canopy than in either gap. Our results highlight that snow accumulation and melt dynamics are controlled by meteorological conditions, soil freezing is controlled by forest structure, and snow properties are controlled by both weather forcing and canopy discontinuity. Overall, observations and simulations suggest that the exceptionally low spring streamflow in the winter of 2020–2120 was mainly driven by low snow accumulation, slow snowmelt, and low precipitation in April and May rather than enhanced percolation through the snowpack and soil freezing.

Modelling, automation, and control are widely used for water resource recovery facility (WRRF) optimization. An influent generator (IG) is a model, aiming to provide the flowrate and pollutant concentration dynamics at the inlet of a WRRF for a range of modelling applications. In this study, a new data-driven IG model is proposed, only using routine data and weather information, and without need for any additional data collection. The model is constructed by an artificial neural network (ANN) and completed with a multivariate regression to generate time series for certain pollutants. The model is able to generate flowrate and quality data (TSS, COD, and nutrients) at different time scales and resolutions (daily or hourly), depending on various user objectives. The model performance is analyzed by a series of statistical criteria. It is shown that the model can generate a very reliable dataset for different model applications.

Abstract A full‐scale biofilm‐enhanced aerated lagoon using fixed submerged media was monitored using automated water quality monitoring stations over the span of one year to quantify its nitrification performance. The system was operating at a high organic loading rate averaging 5.8 g total CBOD 5 /m 2 of media per day (23.9 g total CBOD 5 /m 3 of lagoon per day), a total ammonia nitrogen loading rate averaging 0.9 g NH 4 ‐N/m 2 day (3.7 g NH 4 ‐N/m 3 day), and temperatures ranging from 1.6 to 20.8°C. The system showed an extended seasonal nitrification period compared with a simulated aerated lagoon system of the same dimensions. This extension of complete nitrification with approximately 1 month was observed in the fall despite the decrease of operating temperature down to 4°C. During this maximum nitrification period, substantial denitrification occurred, and the effluent un‐ionized ammonia ratio was reduced. A temporary loss of nitrification was also experienced in relation to an episode of elevated suspended solids concentration. Measured biofilm characteristics, namely the detachment dynamics and the biofilm thickness, were used to explain this temporary nitrification loss. During wintertime, a low nitrate production was still observed, suggesting year‐long retention of nitrifying bacteria in the biofilm. Practitioner points Nitrification in a highly loaded biofilm‐enhanced aerated lagoon is mainly affected by operating temperature. Maximum nitrification is observed during the warmer months and occurs even at high organic loading rates (>5 g CBOD 5 /m 2 day). Compared with a simulated suspended growth system, the biofilm‐enhanced lagoon shows a significantly extended nitrification period. The extension is observed at the end of the summertime maximum nitrification period. Low amounts of nitrate still produced during winter in the biofilm‐enhanced aerated lagoon suggest year‐long retention of autotrophic nitrifying biomass in the biofilm. Nitrification in the biofilm‐enhanced aerated lagoon is negatively impacted by the presence of important quantities of accumulated solids that resuspend when their digestion starts as temperature increases.

A data-driven model was proposed for generating the influent flow and water temperature dynamics including the impact of snowmelt under cold climate conditions. The performance was compared with a phenomenological model.

Abstract. 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). Cold content for a given snowpack thus plays a critical role because it affects both the timing and the rate of snowmelt. Measuring cold content is a labour-intensive task as it requires extracting in situ snow temperature 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–2018. 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 modelled CC was produced and analyzed to highlight the high-resolution temporal variability of CC for the 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 d, with melt events occurring on 42 of those days. Based on snowpit observations, the largest mean CC (−2.62 MJ m−2) 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−2, 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 bulk estimates of snow water equivalent (SWE) (R2=0.64 and percent bias (Pbias) =-17.1 %), snow density (R2=0.71 and Pbias =1.6 %), and cold content (R2=0.93 and Pbias =-3.3 %). A snow density profile derived by utilizing an empirical formulation also provided reasonable estimates of layered cold content (R2=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.

Flood-recession agriculture (FRA) represents a crucial source of livelihood for numerous communities across Africa who reside near expansive floodplains and wetlands. However, it is currently insufficiently monitored. In this study, we present a methodology for mapping FRA harvested areas in the Senegal River Valley that is both reproducible and scalable. Our methodology entails the integration of optical and radar data from Sentinel platforms, conducted through a multitemporal analysis with a seasonal focus, and the application of the Random Forest algorithm. The results, supported by a kappa coefficient of 91.9%, demonstrate the first comprehensive mapping of FRA in the Senegal River valley, conducted between 2019 and 2023. This mapping facilitates the identification of the hydrological factors that influence FRA harvesting. The results of the analyses have demonstrated the importance of interannual variability in the cultivated areas of FRA, which range from 14,000 to 75,000 ha depending on the intensity of the annual flood. The duration and flooded extension are the primary factors that regulate the cropping pattern of FRA over the floodplain. The flood duration must be around 35 days to permit the cultivation, with growth generally starting between 10 and 30 November. In consideration of these findings, we recommend that future water management strategies and rural development initiatives give due consideration to FRA, to enhance the visibility of farmers.

Abstract. In the boreal forest, winter temperatures are projected to increase substantially by 2100, resulting in a reduction in snow cover thickness and duration. These changes are likely to affect hydrological processes such as snowmelt, the soil thermal regime, and snow metamorphism. The exact impact of future changes is difficult to pinpoint in the boreal forest, due to its complex structure, and the fact that snow dynamics under the canopy are very different from those in the gaps. Although the effects of warmer winters on snow-related processes are well documented, their interactions to influence the spring runoff in evergreen forest remain poorly understood. In this observational study, we assess the influence of a low-snow and warm winter on snowmelt dynamics, soil freezing, snowpack properties, and spring streamflow in a humid and discontinuous boreal catchment of eastern Canada (47.29° N, 71.17° W, 850 m ASL). We monitored the soil and snow thermal regimes and sampled physical properties of the snowpack under the canopy and in two forest gaps during an exceptionally low-snow and warm winter, plausibly representative of future winters, and during a winter with conditions close to normal. We observe that snowmelt was earlier but slower, top soil layers were cooler, and gradient metamorphism was enhanced during the low-snow and warm winter. However, we observe that snowmelt duration increased in forest gaps, that soil freezing was enhanced only under the canopy, and that snow permeability increased more strongly under the canopy than in either gap. Overall, we observe that the spring streamflow discharge was significantly reduced in the warmest year due to a slower melt and low precipitation in April and May. Our results, based solely on field observations, highlight the complex effects of warmer winters on snow hydrology in discontinuous boreal forests.

In the alluvial plains of large rivers, annual flooding is essential for numerous ecosystem services, including flood-based agriculture, biodiversity and groundwater recharge. Remote sensing provides increased opportunities to monitor surface water dynamics across large floodplains that are currently poorly captured by local hydrological monitoring and modelling due to data scarcity and the flat, heterogeneous topography. Combining the advances in earth observations with hydrological modelling and extensive in situ fieldwork, this research seeks to improve our understanding of surface water dynamics and associated agricultural practices in the Senegal river floodplain. 2813 mosaics from Landsat, MODIS and Sentinel-2 earth observations are created to map and monitor surface water variations using a site specific MNDWI classification adapted to complex, wetland environments. Validated against ground truth data, the approach is upscaled using cloud computing across this 2250 km 2 floodplain over 1999–2022. Statistical regression models are then developed to estimate flooded and cultivated areas based on upstream flow values since 1950 and analyse trends and exceedance probabilities over time. Results reveal extreme interannual variations in peak flooded areas, ranging from 30,000 ha and 720,000 ha between 1950 and 2022, while annual water modules fluctuate between 210 and 1460 m 3 /s. After 1994, flooded areas show partial recovery, with 95th percentile reaching 89,000 ha during 1994–2022 compared to 37,000 ha in 1972–1993. Flood-based agricultural practices cover between 13,000 ha and 133,000 ha over the same period, highlighting the pronounced variability faced by local rural communities. Occurrence maps and predictive models for annual flooded and cultivated areas based on upstream flows can support early warning tools, helping to prepare for extreme floods and droughts. These outputs are crucial to assess the impact of future climatic and anthropic changes in the region, including planned dams, on the amplitude of annual floods and their associated environmental benefits. • Annual floods in African floodplains support vital flood-based agricultural systems. • Long-term variations monitored through Earth observations & hydrological modelling. • In Senegal River floodplain peak inundated areas reach 300-7200 km 2 over 1950–2022. • Flood recession cropping on 57,000 ha on average since 1994 in this floodplain. • Resulting maps & models support early warning tools and future dam management.

Abstract Accurately modeling the interactions between inland water bodies and the atmosphere in meteorological and climate models is crucial, given the marked differences with surrounding landmasses. Modeling surface heat fluxes remains a challenge because direct observations available for validation are rare, especially at high latitudes. This study presents a detailed evaluation of the Canadian Small Lake Model (CSLM), a one-dimensional mixed-layer dynamic lake model, in reproducing the surface energy budget and the thermal stratification of a subarctic reservoir in eastern Canada. The analysis is supported by multiyear direct observations of turbulent heat fluxes collected on and around the 85-km 2 Romaine-2 hydropower reservoir (50.7°N, 63.2°W) by two flux towers: one operating year-round on the shore and one on a raft during ice-free conditions. The CSLM, which simulates the thermal regime of the water body including ice formation and snow physics, is run in offline mode and forced by local weather observations from 25 June 2018 to 8 June 2021. Comparisons between observations and simulations confirm that CSLM can reasonably reproduce the turbulent heat fluxes and the temperature behavior of the reservoir, despite the one-dimensional nature of the model that cannot account for energy inputs and outputs associated with reservoir operations. The best performance is achieved during the first few months after the ice break-up (mean error = −0.3 and −2.7 W m −2 for latent and sensible heat fluxes, respectively). The model overreacts to strong wind events, leading to subsequent poor estimates of water temperature and eventually to an early freeze-up. The model overestimated the measured annual evaporation corrected for the lack of energy balance closure by 5% and 16% in 2019 and 2020. Significance Statement Freshwater bodies impact the regional climate through energy and water exchanges with the atmosphere. It is challenging to model surface energy fluxes over a northern lake due to the succession of stratification and mixing periods over a year. This study focuses on the interactions between the atmosphere of an irregular shaped northern hydropower reservoir. Direct measurements of turbulent fluxes using an eddy covariance system allowed the model assessment. Turbulent fluxes were successfully predicted during the open water period. Comparison between observed and modeled time series showed a good agreement; however, the model overreacted to high wind episodes. Biases mostly occur during freeze-up and breakup, stressing the importance of a good representation of the ice cover processes.

At high latitudes, lake-atmosphere interactions are disrupted for several months of the year by the presence of an ice cover. By isolating the water column from the atmosphere, ice, typically topped by snow, drastically alters albedo, surface roughness, and heat exchanges relative to the open water period, with major climatic, ecological, and hydrological implications. Lake models used to simulate the appearance and disappearance of the ice cover have rarely been validated with detailed in situ observations of snow and ice. In this study, we investigate the ability of the physically-based 1D Canadian Small Lake Model (CSLM) to simulate the freeze-up, ice-cover growth, and breakup of a small boreal lake. The model, driven offline by local weather observations, is run on Lake Piché, 0.15 km 2 and 4 m deep (47.32°N; 71.15°W) from 25 October 2019 to 20 July 2021, and compared to observations of the temperature profile and ice and snow cover properties. Our results show that the CSLM is able to reproduce the total ice thickness (average error of 15 cm) but not the ice type-specific thickness, underestimating clear ice and overestimating snow ice. CSLM manages to reproduce snow depth (errors less than 10 cm). However, it has an average cold bias of 2°C and an underestimation of average snow density of 34 kg m −3 . Observed and model freeze-up and break-up dates are very similar, as the model is able to predict the longevity of the ice cover to within 2 weeks. CSLM successfully reproduces seasonal stratification, the mixed layer depth, and surface water temperatures, while it shows discrepancies in simulating bottom waters especially during the open water period.

Abstract. 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 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 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 observations, the largest mean CC (−2.62 MJ m−2) 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−2, 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) (R2 = 0.64 and percent bias (Pbias) = −17.1 %), bulk snow density (R2 = 0.71 and Pbias = 1.6 %), and bulk cold content (R2 = 0.90 and Pbias = −2.0 %). A snow density profile derived by utilizing an empirical formulation also provided reasonable estimates of cold content (R2 = 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.

Abstract. 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 leads 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 favored 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 27 of the 32 observed melt–freeze layers induced by rain-on-snow events to be formed by the model, instead of only 18 with the original canopy module. This new model development also delays and reduces the snowpack runoff. Indeed, it triggers the unloading of dense unloaded snow layers with small rounded grains, which in turn produces fine–over–coarse transitions that limit percolation and favor refreezing. Our results show 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 properties measurements to improve hydrological models in forested snow–covered regions.