Benjamin Bouchard
email: benjamin.bouchard.1@ulaval.ca

Daniel F Nadeau

Florent Domine

François Anctil

Tobias Jonas

Étienne Tremblay

How does a warm and low-snow winter impact the snow cover dynamics in a humid and discontinuous boreal forest? An observational study in eastern Canada

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The boreal forest is one of the most extensive biomes on Earth. It is projected to warm by up to 5°C by 2100, with the largest increases occurring in winter [START_REF] Scheffer | Thresholds for boreal biome transitions[END_REF][START_REF] Zheng | Record-high CO2 emissions from boreal fires in 2021[END_REF]Price et al., 2013;IPCC, 2022). Warmer winters will result in less solid precipitation and in a thinner, shorter-lived snow cover [START_REF] Laternser | Long-term snow climate trends of the Swiss Alps (1931-99)[END_REF][START_REF] Hamlet | Effects of temperature and precipitation variability on snowpack trends in the Western United States[END_REF]. Spring melt will occur earlier in the season, but at a slower rate because less radiative energy is available at that time [START_REF] López-Moreno | Response of snow processes to climate change: spatial variability in a small basin in the Spanish Pyrenees[END_REF][START_REF] Musselman | Slower snowmelt in a warmer world[END_REF]. Together, these changes are expected to reduce peak spring streamflow and runoff volumes [START_REF] Furey | A stochastic conceptual modeling approach for examining the effects of climate change on streamflows in mountain basins[END_REF]Berghuijs et al., 2014;[START_REF] Luce | Declining annual streamflow distributions in the Pacific Northwest United States, 1948-2006[END_REF][START_REF] Barnhart | Snowmelt rate dictates streamflow[END_REF]. To date, most studies focusing on hydrological changes associated with warmer winters in the boreal forest have looked at rather dry regions of the biome, in western Canada and Europe [START_REF] Ireson | The changing water cycle: the Boreal Plains ecozone of Western Canada[END_REF][START_REF] Shook | Changes in the hydrological character of rainfall on the Canadian prairies[END_REF][START_REF] Xu | The Effect of Climate Change on River Flow and Snow Cover in the NOPEX Area Simulated by a Simple Water Balance Model[END_REF]. The limited number of studies on the humid boreal forest of eastern Canada point towards an increase in mean annual streamflow, but a decrease in snow accumulation and an earlier spring freshet [START_REF] Guay | A global portrait of hydrological changes at the 2050 horizon for the province of Québec[END_REF]Ouranos and MELCCFP, 2022;[START_REF] Valencia Giraldo | Assessment of the Potential Hydrological Impacts of Climate Change in Quebec-Canada, a Refined Neutral Approach[END_REF].

The magnitude of the expected hydrological changes in response to warmer winters depends on a number of processes that occur in the soil-snow-forest-atmosphere continuum. One such process is soil freezing. In forests, the spatial pattern of soil freezing is difficult to determine because snow depth is highly variable and freezing depth and snowpack thickness are inversely proportional [START_REF] Goodrich | The influence of snow cover on the ground thermal regime[END_REF][START_REF] Zhang | Influence of the seasonal snow cover on the ground thermal regime: An overview[END_REF][START_REF] Hardy | Snow depth manipulation and its influence on soil frost and water dynamics in a northern hardwood forest[END_REF]. Therefore, it is challenging to understand and predict infiltration versus surface runoff in forested environments during winter. Moreover, the heterogeneity of forested soil freezing is increased by the abundance of macropores [START_REF] Beven | Macropores and water flow in soils revisited[END_REF]. On the one hand, macropores initially rich in air create preferential pathways for infiltration [START_REF] Granger | Snowmelt infiltration to frozen Prairie soils[END_REF][START_REF] Stähli | Preferential water flow in a frozen soil -A two-domain model approach[END_REF][START_REF] Fuss | Hydrologic flowpaths during snowmelt in forested headwater catchments under differing winter climatic and soil frost regimes[END_REF]. On the other hand, macropores that are rich in ice, for example following a refreeze [START_REF] Kane | Water movement into seasonally frozen soils[END_REF][START_REF] Demand | Influences of macropores on infiltration into seasonally frozen soil[END_REF], significantly reduce further infiltration [START_REF] Mohammed | Snowmelt infiltration and macropore flow in frozen soils: Overview, knowledge gaps, and a conceptual framework[END_REF][START_REF] Watanabe | Effect of macropores on soil freezing and thawing with infiltration[END_REF]. Previous studies have reported that both canopy interception, which limits snowpack thickening, and meteorological conditions such as rain-on-snow events followed by a cold spell, which favor water refreezing within macropores, promote runoff over infiltration [START_REF] Proulx | Classification of meteorological conditions to assess the potential for concrete frost formation in boreal forest floors[END_REF]Jones and Pomeroy, 17-19 May 2001;[START_REF] Stadler | Measurement of frost-induced snowmelt runoff in a forest soil[END_REF][START_REF] Shanley | The effect of frozen soil on snowmelt runoff at Sleepers River, Vermont[END_REF]. With the projected warmer climate leading to more frequent and intense winter rainfalls (Il Jeong and Sushama, 2018), it is therefore important to deepen our understanding of soil freezing dynamics in discontinuous boreal forests.

Forest structure affects not only the dynamics of soil freezing, but also the physical properties of the snowpack. As canopy interception limits snow accumulation under trees [START_REF] Pomeroy | Coupled modelling of forest snow interception and sublimation[END_REF][START_REF] Mazzotti | Revisiting snow cover variability and canopy structure within forest stands: Insights from airborne LiDAR data[END_REF][START_REF] Sun | Evaluating the functionality and streamflow impacts of explicitly modelling forest-snow interactions and canopy gaps in a distributed hydrologic model[END_REF], stronger vertical temperature gradients (𝜕T/𝜕z) favor kinetic snow grain growth under the canopy more than in gaps [START_REF] Albert | Snowpack stratigraphy evolution at forested and open sites[END_REF][START_REF] Molotch | Measuring spatiotemporal variation in snow optical grain size under a subalpine forest canopy using contact spectroscopy[END_REF][START_REF] Bouchard | Comparison of snowpack structure in gaps and under the canopy in a humid boreal forest[END_REF]. [START_REF] Bouchard | Comparison of snowpack structure in gaps and under the canopy in a humid boreal forest[END_REF] showed that the phenomenon results in lower specific surface area (SSA) and greater snowpack permeability under the canopy. With warmer winters, one may expect snow surface temperature to increase, but snow thickness to decrease. [START_REF] Domine | Interactions between snow metamorphism and climate: Physical and chemical aspects[END_REF] have suggested a decrease in permeability with climate warming due to a warmer snow surface, lower 𝜕T/𝜕z, and slower grain growth. The authors also 3 noted that an increase in the frequency of melting events would result in more melt-freeze crusts and low-permeability ice layers [START_REF] Albert | Ice layer and surface crust permeability in a seasonal snow pack[END_REF]. In the boreal forest, where snow surface temperature and snow accumulation are highly variable in space [START_REF] Malle | Shading by trees and fractional snow cover control the subcanopy radiation budget[END_REF]Parajuli et al., 2020;Mazzotti et al., 2023) and where the canopy structure influences crust formation in the snowpack [START_REF] Bründl | Routing of canopy drip in the snowpack below a spruce crown[END_REF][START_REF] Teich | Effects of bark beetle attacks on forest snowpack and avalanche formation -Implications for protection forest management[END_REF][START_REF] Bouchard | Comparison of snowpack structure in gaps and under the canopy in a humid boreal forest[END_REF], the impacts of warmer winters on snow properties may be hard to predict. Therefore, dedicated studies on the influence of warmer winters on snowpack physical properties in discontinuous boreal forests are needed.

Based on the research gaps identified above, we defined two objectives in this study. The first one is to quantify the effect of a low-snow and warm winter on snowmelt, soil freezing and snowpack physical properties in a humid and discontinuous boreal forest. The second one is to evaluate how these factors interact together to modulate spring runoff. To assess this, we compared snow melt, snow physical properties, soil freezing, and spring runoff at a small catchment in the humid boreal forest of eastern Canada, for two consecutive winters. One winter was exceptionally warm and dry, while the other was slightly colder, with precipitation amounts similar to the standard climatology of the study region. These contrasted conditions represent an ideal comparison to investigate some expected effects of climate change. Extensive snow monitoring and pit measurements were performed to contribute to the objectives.

Section 2 presents the study site, the instrumentation, and the observed and estimated physical variables. In Sect. 3, we present the climatology of each winter and the differences in melt rate, snowpack and ground thermal regime, and the evolution of snow cover physical properties within medium-size and small forest gaps and under the canopy. The resulting spring streamflow of the catchment is also presented for both winters in Sect. 3. Finally, we compare our results with existing literature, present the limitations of the study, and discuss the potential climatic, hydrological, and ecological implications of warmer winters in a discontinuous humid boreal forest in Sect. 4.

Methods

Study site

The study took place in the Montmorency Forest (MF) during winters 2020-2021 (W20-21) and 2021-2022 (W21-22).

Measurements began on 15 October and ended on 15 June of the following year. MF is located in the province of Québec, in eastern Canada, at the southern edge of the boreal forest (47.29° N, 71.17° W; Fig. 1a). Three monitoring stations were installed in a 1.1 km 2 forested subcatchment (7A) of the Bassin expérimental du Ruisseau des Eaux-Volées (BEREV) located at MF (Fig. 1b). The stations were installed in a medium gap, a small gap, and under the canopy (Fig. 1c). Table 1 presents the gap dimensions and the sky-view fraction (SVF) at each site. The SVF was evaluated using an adaptive thresholding algorithm [START_REF] Jonas | HPEval: A canopy shortwave radiation transmission model using highresolution hemispherical images[END_REF] on hemispherical photographs taken with a Sigma 8mm f3.5 EX DG Circular Fisheye lens in the fall of 2022 (Fig. 1d-e-f). Stations were located on a 12°, northeast-facing slope at 846 m above sea level (ASL), within a 10 m tall stand of balsam fir (Abies balsamea) mixed with white birch (Betula papyrifera) and white spruce (Picea glauca) trees. The catchment is a regeneration from a major logging operation that took place in 1993 [START_REF] Guillemette | Rainfall generated stormflow response to clearcutting a boreal forest: peak flow comparison with 50 world-wide basin studies[END_REF]. The soil is a sandy loam topped by »7 cm of litter. The stations were located in the vicinity of a 15 m flux tower measuring shortwave (0.3-2.8 µm) and longwave (4.5-42 µm), upwelling and downwelling radiation (CNR4;

Kipp & Zonen). A V-notch streamflow gauge and a bubble flowmeter, maintained by the provincial government and in operation since 1967, are located at the outlet of the catchment (DEH station 051004, https://www.cehq.gouv.qc.ca/hydrometrie/historique_donnees/). Based on LiDAR imagery, we estimate that 75 % of the catchment is canopy-covered, with the remainder being gaps of small and medium size. Four kilometers northeast of BEREV-7A, there is a 0.01 km 2 open area at 664 m ASL, the NEIGE site [START_REF] Pierre | Evaluation of catch efficiency transfer functions for unshielded and single-alter-shielded solid precipitation measurements[END_REF], which hosts a federal weather station (ECCC station 7042395). We used data from the federal station to position the winters of 20-21 and 21-22 in relation to the MF climatology. The NEIGE site also hosts a CS725 instrument (Campbell Scientific) that monitors the snow water equivalent (SWE) at a 6 h timestep, using differential gamma-ray absorption [START_REF] Choquette | CS725, an accurate sensor for the snow water equivalent and soil moisture measurements[END_REF]) and an ultrasonic snow height sensor (Judd Communication) working on an hourly timestep.

Measurements of physical variables

In this study, the monitoring stations tracked the thermal regime of the snowpack and of the top 20 cm of soil, the soil volumetric water content (VWC) and the effective thermal conductivity of the snow (ks) at two levels in the snowpack.

Measurements from the stations were complemented with snow pit observations of snow density (𝜌 ! ), temperature, SSA, and a visual identification of the snowpack stratigraphy conducted periodically during both winters.

Monitoring Stations

Each station includes a snow height sensor, either the Judd model or the SR50a from Campbell Scientific, mounted 3 m above the ground, and is equipped with a vertical array of PT1000 thermistors (Schneider Electric). Temperature was measured at depths of 20, 10 and 5 cm in the soil, at the ground surface and then every 15 cm in the snowpack until a maximum height of 180 cm. Thermistors were inserted into white-painted aluminum tubes held in place by a vertical UHMW plastic rod. Snow surface temperature was measured at each station with an SI-111 infrared radiometer (Apogee Instruments) mounted at the same level as the snow height sensor. Each station was also equipped with an HMP60 (Campbell Scientific) sensor that measured air temperature and relative humidity at a height of 3 m above the ground. A CS655 reflectometer (Campbell Scientific) measured the VWC at a depth of 15 cm into the soil. Variations in 15 cm VWC were used as a proxy for infiltration.

A CR10X data logger (Campbell Scientific) recorded point measurements every hour. Also, time-lapse cameras taking hourly photos were used to visually identify the precipitation phase and the burying of the thermistors. An example of a monitoring station, the one under the canopy in this case, is presented in Fig. 1g.

The stations were also equipped with two heated TP08 needle probes (NP; Hukseflux) to measure ks following Morin et al.

(2010) and Domine et al. (2015). The lowest NP was installed 10 cm above the ground at each station, whereas the highest NP was installed 65 cm above the ground under the canopy, and 80 cm in the forest gaps to account for the greater snow height in these environments. To prevent snow melting around the heated needle during measurements, these were taken only every other day between 5:00 and 6:00 AM, only when the snow temperature was below -2.5°C. During each measurement, the NP was heated for 150 s with a power of 0.45 W m -1 . The heating curve was recorded on a CR1000 data logger (Campbell Scientific). The algorithm developed to automatically assess ks from the heating curves is derived from Domine et al. (2015) https://doi.org/10. and is described in the Supplementary Material (Fig. S1). The error associated with snow thermal conductivity measurements with these static needle probes is estimated to be 21 % (Domine et al., 2015).

Data gap filling

Small data gaps (1 to 5 h) in the snow height and temperature time series were filled by linear interpolation. We used data from the other stations to fill the longer data gaps, as described in the Supplementary Material (Fig. S2 -S4). This approach was applied to snow height in the medium forest gap (37 days in 2020) and under the canopy (77 days in 2022), to snow surface temperature in the small forest gap (winter 2021-22), and to air and snow surface temperatures under the canopy (16 days in 2022).

Snow pit measurements

Four snow pits were dug in medium gaps, small gaps and under the canopy each winter, for a total of 24 snow pits over the two winters. Snow pits in gaps and under the canopy were all dug at sites with similar conditions (gap size, SVF) to the monitoring stations, all within 150 m of the stations. Table 2 lists the date of each snow pit measurement.

For each snow pit survey, we measured 𝜌 ! with a Snow-Hydro 100 cm 3 box with a ±10% accuracy [START_REF] Conger | Comparison of density cutters for snow profile observations[END_REF]. Snow density was measured every 3 to 5 cm in the vertical direction. In the presence of ice columns, measurements were taken adjacently. We measured snow temperature using a Greinsinger PT-1000 probe (resolution: 0.1 °C) every 5 cm in the topmost 40 cm of snow and every 10 cm for the lower snow layers. The SSA, being a quantitative indicator of snow metamorphism [START_REF] Domine | Interactions between snow metamorphism and climate: Physical and chemical aspects[END_REF], was measured using the DUFISSS instrument [START_REF] Gallet | Measurement of the specific surface area of snow using infrared reflectance in an integrating sphere at 1310 and 1550 nm[END_REF]. DUFISSS uses infrared reflectance of snow samples with an integrating sphere at 1310 nm to estimate the SSA with an accuracy of approximately 12% [START_REF] Gallet | Measurement of the specific surface area of snow using infrared reflectance in an integrating sphere at 1310 and 1550 nm[END_REF]. Measurements were taken vertically every 1 to 5 cm depending on the stratigraphy.

Estimated variables

Snow water equivalent and snowpack cold content

We estimated the SWE (m) from snow density profiles using:

𝑆𝑊𝐸 = "# ! $$$ # " , (1)
where ℎ and 𝜌 ! ) are respectively the height (m) and the average density (kg m -3 ) of the snowpack, and 𝜌 % is the density of liquid water (1000 kg m -3 ).

The snowpack cold content (QCC, in J m -2 ) is the quantity of energy needed to bring the snowpack to its melting point. It is defined as:

𝑄 && = -𝑐 '() ∑ ℎ ' 𝜌 !,+ .... (𝑇 !,+ .... -𝑇 , ) - './ , ( 2
)
where 𝑐 '() is the heat capacity of ice (2108 J kg -1 K -1 ), ℎ ' and 𝜌 !,+ .... are the same as above but for a given snow layer 𝑖. 𝑇 !,+ .... is the average temperature between the bottom and top boundaries of each 𝑖 layer (K) as monitored at the stations and 𝑇 , is the melting point of ice (273.15 K). ℎ ' corresponds to the vertical distance between each temperature measurements.

Since snow density profiles were point measurements, we interpolated 𝜌 !,+ .... linearly between each snow pit date at an hourly timestep to match temperature monitoring. Snow density was assumed constant from the beginning of the snow season until the date of the first snow pit (23 days in W20-21 and 59 days in W21-22) and then from the last snow pit date until the meltout (39 days in W20-21 and 31 days in W21-22). To validate this rough hypothesis, we measured four density profiles in W20-21 and five in W21-22 at NEIGE site next to the SWE sensor and compared the resulting estimated SWE time series with direct observations, as shown in the Supplementary Material (Fig. S5).

Snowpack net shortwave and longwave radiation

We used the HPEval model [START_REF] Jonas | HPEval: A canopy shortwave radiation transmission model using highresolution hemispherical images[END_REF] to estimate the downwelling shortwave radiation below the canopy (𝑆𝑊𝑅 ↓,1( ) from the hemispherical photographs taken at the monitoring stations and the incoming shortwave radiation measured above the canopy at the flux tower, some 10 m away. HPEval combines hemispherical imagery of the canopy and radiation transfer modeling to estimate subcanopy shortwave radiation at very high spatial and temporal resolution. Reflected shortwave radiation by the snowpack (𝑆𝑊𝑅 ↑,1( ) was estimated using five arbitrary albedo classes defined from the work from [START_REF] Hardy | Incorporating effects of forest litter in a snow process model[END_REF] and [START_REF] Melloh | Spectral albedo/reflectance of littered forest snow during the melt season[END_REF]. The five classes are listed in Table 3.

We manually assigned a daily albedo class to the snow surface by visually inspecting timelapse photographs from each station.

A sample representative photo of each albedo class is presented in Figure S6. The net shortwave radiation below canopy (𝑆𝑊𝑅 -)3,1( ) is the difference between 𝑆𝑊𝑅 ↓,1( and 𝑆𝑊𝑅 ↑,1( . We estimated the downwelling longwave radiation below the canopy (𝐿𝑊𝑅 ↓,1( ; W m -2 ) by using:

𝐿𝑊𝑅 ↓,1( = 𝑆𝑉𝐹 × 𝐿𝑊𝑅 ↓,4( + (1 -𝑆𝑉𝐹) 𝜀 (4-𝜎 𝑇 (4- 5 , (3)
where 𝐿𝑊𝑅 ↓,4( is the downwelling longwave radiation measured above the canopy at the nearby flux tower (W m -2 ), 𝜎 is the Stefan-Boltzmann constant (5.67 × 10 -8 W m -2 K -4 ), 𝜀 (4-is the canopy emissivity, set to 0.98 (Pomeroy et al., 2009), and Tcan is the canopy temperature (K). We assume that the canopy temperature can be approximated by the air temperature measured at our monitoring stations.

The upwelling longwave radiation below the canopy (𝐿𝑊𝑅 ↑,1( ; W m -2 ) is determined with:

𝐿𝑊𝑅 ↑,1( = 𝜀 ! 𝜎 𝑇 !678 5 , ( 4
)
where 𝜀 ! is the emissivity of the snow surface, set to 0.99 (Sicart et al., 2006), and Tsurf is the snow surface temperature (K) measured at each station. The net longwave radiation below canopy (𝐿𝑊𝑅 -)3,1( ) is the difference between 𝐿𝑊𝑅 ↓,1( and 𝐿𝑊𝑅 ↑,1( . The net total radiation below canopy (𝑅 -)3,1( ) is the sum of 𝑆𝑊𝑅 -)3,1( and 𝐿𝑊𝑅 -)3,1( .

Ground heat flux

The ground heat flux was assumed equivalent to the snow heat flux within the lower 15 cm of the snowpack (Lackner et al., 2022), and calculated using Fourier's law (Eq. 5).

𝐹 = -𝑘 ! (: #$ ;: % ) =" , (5)
where F (W m -2 ) is the heat flux of the lower 15 cm of the snowpack, ks (W m -1 K -1 ) is the effective thermal conductivity of this basal snow layer, taken as the NP measurement at a 10 cm height, T15 and T0 are the temperatures measured 15 cm above the ground and at the ground surface (K), respectively, dh is the thickness of the bottom layer (15 cm). To avoid melting of the ice matrix, NPs were not heated when the snow was warmer than -2.5°C. Because of this constraint, ks at 10 cm from each station was only measured at the beginning of the winter on both years when the bottom of the snowpack was colder than -2.5°C. This resulted in more than 80% ks measurements that were missing. In order to have complete time series of ks, we used Eq. 18 from Fourteau et al. (2021) with our pit density measurements to estimate ks. We validated the use of this empirical equation by comparing it to our observations of pit density and monitoring of ks at 10 and 65 or 80 cm at all stations. We used ks that was measured the closest in time to the snow pit surveys. Before comparing our observations with the equation, correction factors between 1.1 and 1.3, depending on snow temperature, were applied to ks measurements to account for the systematic underestimation of the fixed NP approach (Fourteau et al., 2022). Doing this, we found a correlation coefficient of 0.70 and a bias of 10.2% between our observations and the equation from Fourteau et al. (2021). Comparison between the corrected measurements and the equation is detailed in the Supplementary Material (Fig. S7).

Vertical temperature gradient

The magnitude of the mean snowpack vertical temperature gradient (|𝜕T/ 𝜕z|; in °C m -1 ) is expressed as follow:

|𝜕𝑇/ 𝜕𝑧| = @ >: !&'( ;: % ? " ! @ , (6)
where (Tsurf -T0) is the temperature difference between the snow surface and the soil-snow interface (in °C) and the hs is the snow height (in m).

Snow permeability

Snow permeability (Ks; in m 2 ) indicates the ease with which a fluid subjected to a pressure gradient flows through a porous medium. Ks can be estimated from 𝜌 ! and the optical grain radius (rg; in m) according to the following equation (Calonne et al., 2012):

𝐾 ! = 3𝑟 @ A 𝑒 ;B.B/D# ! . ( 7
)
We assume that each snow layer is a collection of independent ice spheres (𝜌 '() = 917 kg m -3 ) with the same SSA as the snow of interest [START_REF] Grenfell | Representation of a nonspherical ice particle by a collection of independent spheres for scattering and absorption of radiation[END_REF]. Therefore, rg is estimated as follow:

𝑟 @ = D # )*+ EEF . ( 8
)
3 Results

Climatic conditions

Figure 2 shows the monthly sum of precipitation and the mean of the monthly air temperature from October to May, as recorded at the NEIGE site in W20-21 and W21-22 compared with the average over 1982-2022. W20-21 was the driest winter of the 1982-2022 period, with 199 mm recorded from January to April (JFMA), including 167 mm of solid precipitation (Table 4).

This corresponds to a precipitation anomaly of -224 mm. In comparison, the JFMA anomaly of precipitation in W21-22 was only +5 mm. December-January-February (DJF) temperature in W20-21 was also 2.8°C warmer than the 1982-2022 average.

As for DJF, W20-21 was the fourth warmest year of the last 40 years. In comparison, DJF temperature was 1.8°C colder in W21-22 than the 1982-2022 average and ranks as the sixth coldest for that period. The exceptionally dry and warm winter of 4). Two heavy rain-on-snow (ROS) events were observed in December 2020, three smaller ones in November and December 2021, and several others in March, April, and May of both years. As the data show, the winter 2020-21 received a much lower than average solid precipitation and was significantly warmer resulting in a thin snowpack. We will refer to it as our low-snow and warm winter. The precipitation and temperature anomaly in the winter of 2021-22 is much weaker. Therefore, W21-22 will correspond to the reference winter for the rest of the analysis.

Snow accumulation and melt dynamics

Figure 3 shows the evolution of the snowpack thermal regime in gaps and under the canopy, along with the snowpack cold 255 content (QCC). Due to a lower snowfall, the snowpack in W20-21 was on average 46% thinner than in the reference winter.

Since the air was warmer in W20-21, the temperature of the upper snow layers was also higher. In contrast, because there was less snow on the ground in W20-21, heat transfer through the snow was facilitated, resulting in the base of the snowpack being colder in W21-22. Due to canopy interception, less snow accumulated under the trees than in the gaps in both years.

Interestingly, topmost snow layers seem to be colder under the canopy than inside gaps. This is in contradiction with Figure 4 260 that shows snow surface temperature being higher under the canopy than inside gaps at night, despite similar air temperature at all three stations. This discrepancy is further discussed in Sect. 4.2. Due to the thinner snowpack and warmer air temperature, the QCC was on average 36% lower in W20-21 than in W21-22. 270 Interestingly, the QCC peaked in late February in both winters. However, since less energy was required to warm the snowpack to 0°C, snowmelt started on average 23 days earlier in the low-snow and warm winter than in the reference winter (10 Avril 2021 vs 3 May 2022), which is a substantial difference. Since the snowpack was thicker in the gaps than under the canopy, the QCC was slightly less under the trees than in both gaps.

As snowmelt started more than three weeks earlier in W20-21 than in the reference winter, less radiative energy was available 275 to contribute to snowmelt (Fig. 5). Net shortwave (𝑆𝑊𝑅 -)3,1( ) and longwave (𝐿𝑊𝑅 -)3,1( ) radiation below canopy were both lower in the low-snow and warm year. As expected, we observed a decrease in 𝑆𝑊𝑅 -)3,1( and an increase in 𝐿𝑊𝑅 -)3,1( as the sky-view-fraction decreased regardless of the year. Figure 6 shows the evolution of SWE from the beginning of snowmelt in both years. In the medium gap, small gap and under the canopy, the SWE at the beginning of snowmelt in the low-snow year was 56%, 59% and 76% lower, respectively, than in the reference year. In addition to a thinner snowpack, the melt rate was also significantly smaller during the low-snow winter. 285

Overall, from W20-21 to W21-22, the duration of the melt period slightly increased from 26 to 31 days and from 29 to 37 days in the medium and small gaps, respectively, whereas it did not change under the canopy.

Ground thermal regime and water content

Figure 7 shows the ground heat flux (GHF) and the soil temperature at several depths in both years. Overall, the GHF in DJF was on average 50% higher in W20-21 than in the reference winter. The largest difference was observed under the canopy, where the GHF was larger than in gaps in W20-21 and in W21-22.

At any given depth, the soil was cooler in W20-21 than in W21-22 in DJF (Fig. 7; Table 5), even though it was a warmer season. In both winters, the top few centimeters of soil below the canopy dropped below freezing, while all soil layers in the gaps remained above or at 0°C. We also observed a frozen soil-snow interface during snow pit surveys under the canopy, but never in the gaps. In the low-snow and warm winter, the freezing depth reached 10 cm, compared to only 5 cm in the reference winter, even though it was much colder. However, the soil-snow interface was below 0°C for similar durations (103 days vs 102 days). Figure 8 shows that soil volumetric liquid water content (VWC) at 15 cm depth was higher under the canopy than in gaps for both years, except in spring and after two heavy ROS that occurred early in W20-21. Soil profile characterization that was performed in a small gap and under the canopy during the summer of 2020 showed different porosity in the topmost 30 cm of soil which could explain the differences in VWC over a short distance at our study site. On 10 April 2021 (Fig. 8a), the increase in VWC was sharp and sudden in both gaps, while it was more gradual under the canopy. A data acquisition error occurred at the canopy station in April and May 2022 (Fig. 8b), so we cannot compare the increase in 15 cm VWC under the canopy at the onset of snowmelt in W20-21 and W21-22. In the medium and small gaps, we observed a similar behavior of the 15 cm VWC at the beginning of snowmelt in 2022.

Vertical temperature gradient and snow properties

Figure 9 shows the magnitude of the vertical snow temperature gradient (|𝜕T/𝜕z|) between the snow surface and the ground.

|𝜕T/𝜕z| decreased throughout the season in both winters. From November to January, the average |𝜕T/𝜕z| was similar at each site and year. In February and March of the low-snow and warm winter, |𝜕T/𝜕z| remained within or above the transition zone from equilibrium to kinetic crystal growth (Colbeck, 1983). In contrast, in February of the reference winter, the decrease of |𝜕T/𝜕z| was more intense and we observed a drop of |𝜕T/𝜕z| below 10°C m -1 in gaps. |𝜕T/𝜕z| was also higher under the canopy than in gaps. At all sites, there was a greater proportion of faceted crystals (FC) and depth hoar (DH) in the snowpack during the low-snow year than during the reference year. In contrast, we observed fewer rounded grains (RG) in W20-21 than in W21-22. In both years, FC and DH layers were proportionally thicker under the canopy than in the gaps. A thick layer of melt-freeze polycrystals (MFpc) was observed at the base of the snowpack in W20-21 as a result of the December 2020 ROS. Also, in W20-21, we observed ice layers (IL) and melt-freeze crusts (MFcr) in the gaps, but not under the canopy. In W21-22, a thin

MFpc was observed at the base of the snowpack. Three to four IL and MFcr were also observed within the gap and canopy snowpacks. Snow density was lower in W20-21 than in the reference winter (257 kg m -3 versus 276 kg m -3 on average). In general, the density profiles showed values increasing with depth and were lower under the canopy than in gaps in both winters. A Kruskal-Wallis test performed on density measurements from all sites showed that the density difference between both winters was significant (p-value < 0.05).

The vertical profile of SSA followed the expected general shape, with higher values near the surface (fresh snow) and lower values deeper (FC and DH). In the bottom 25% of the snowpack, there was no significant difference in SSA between both winters (p-value < 0.05). This may be due to the ice and melt-freeze layers SSA measurement technique that was improved from W20-21 to W21-22. Removing the ice and melt-freeze layers from the analysis leads to a significantly greater SSA in W21-22, coherent with a lower temperature gradient. At a normalized height of 0.25 to 0.75, where the transition from RG to FC occurred, the SSA was 19% lower during the low-snow winter than in the reference winter, which was statistically significant (p-value < 0.05). As a result of a lower density and SSA, the snow permeability was 57% higher in the warm and low-snow winter within the height range 0.25 to 0.75. This difference between both years was also statistically significant, as demonstrated by a Kruskal-Wallis test (p-value < 0.05). Note that in both years, the SSA and the snow permeability was lower and higher, respectively, under the canopy than in gaps.

Spring streamflow

Figure 11 shows the evolution of air temperature and precipitation from the NEIGE site, the decrease in SWE estimated at the monitoring stations and the hydrographs from measurements of spring streamflow at the outlet of the BEREV-7A catchment in April and May of both years. In 2021, air temperatures became positive in early April, allowing a decline in SWE. In 2022, this occurred 20 days later. The daily melt rate was much smaller in the first winter than in the second, as already shown in Fig. 6. The April and May precipitation was also lower in 2021 compared to the following year (122 mm vs 253 mm).

As shown in Fig. 11e-f, spring discharge began to increase on April 7 during the low-snow and warm winter, and began to rise on 1 May, 24 days later, during the reference winter. The spring freshet was earlier than the median spring hydrograph in W20-21, and slightly later in W21-22. The peak discharge in the low-snow winter was similar to the peak of the median hydrograph but was a factor of 5 lower than in W21-22. In both winters, the peak discharge was measured during a day with snowmelt and heavy precipitation (Fig. 11c-d). Low snowmelt and precipitation water inputs in the first winter resulted in significantly lower spring runoff volumes than in the reference winter. In April and May 2021, 2.1 ´ 10 5 m 3 flowed out of the catchment compared to 5.7 ´ 10 5 m 3 in the following year. The spring runoff volume from the low-snow and warm winter was the lowest volume observed at the outlet of BEREV-7A for April and May since discharge monitoring began in 1968. In 2022, it was the sixth highest.

Low-snow and warm winter

In eastern Canada, as in high-latitude and high-altitude regions, less solid precipitation and higher temperatures are expected in winter with climate change, leading to a reduction in snow cover extent [START_REF] Guay | A global portrait of hydrological changes at the 2050 horizon for the province of Québec[END_REF]Ouranos and MELCCFP, 2022;Pepin et al., 2015;Kunkel et al., 2016). Winter 2020-21 at MF received 211 mm less solid precipitation and was warmer by 4.6°C in DJF, and snow melt-out occurred 34 days earlier than in 2021-22, which is more representative of the prevailing local climatic conditions (Fig. 2; Table 4). Indeed, W20-21 was exceptional in that it had both the lowest snowfall in the last 40 years and was one of the warmest winters in that period (ranked fourth). In contrast, W21-22 was a relatively cold winter with an anomaly of -1.8°C compared to the 1982-2022 period, with solid precipitation close to the climate normals. Another feature of global warming in eastern Canada is an increase in the occurrence and intensity of ROS events [START_REF] Vincent | Changes in daily and extreme temperature and precipitation indices for Canada over the twentieth century[END_REF]Il Jeong and Sushama, 2018). Although there were fewer ROS events in W20-21 than in the following year before March (2 versus 3), they were by far more intense, with 44 and 106 mm of liquid precipitation each. These two events reduced the snow cover thickness and triggered winter snowmelt as shown in Fig. 3a-c

Earlier and slower melt

During the low-snow and warm winter, snowpack QCC was lower than in the reference winter. This is partly explained by a thinner snowpack with lower density (Fig. 3 and Fig. 10b-f-j). Warmer upper snow layers in W20-21 also contributed to the lower QCC. However, the lower snow layers were cooler in W20-21, which reduced the difference in QCC between the two years. Our observations are consistent with simulations by [START_REF] Jennings | Snowfall fraction, cold content, and energy balance changes drive differential response to simulated warming in an alpine and subalpine snowpack[END_REF] for alpine and subalpine sites in the western US, where QCC is expected to decrease with increasing winter temperatures. Interestingly, the subcanopy snowpack had a larger QCC than in the forest gaps, despite being thinner and lighter. This would imply an enhanced outgoing heat flux under the canopy that would considerably cool down the snowpack compared to both gaps. However, as noted in section 3.2, this explanation would contradict the warmer snow surface under the canopy in both years (Figure 4). This behavior is expected as radiative cooling in forest is more likely to take place at locations of higher sky-view fraction, such as gaps. This suggests a cold bias in the monitored temperature profile time series under the canopy which lead to an overestimation of the cold content. Timelapse images taken at the canopy station indicate that the snowpack seemed to grow thicker under the ultrasonic sensor than at the thermistor array in both years, which could explain the discrepancy between Figure 3e-f and Figure 4.

Since the snowpack had a lower QCC in W20-21 in both the gaps and below the canopy, less energy was required to bring it to the melting point so snowmelt started earlier, consistent with what is expected in warmer winters [START_REF] Barnett | Potential impacts of a warming climate on water availability in snowdominated regions[END_REF][START_REF] López-Moreno | Response of snow processes to climate change: spatial variability in a small basin in the Spanish Pyrenees[END_REF][START_REF] Guay | A global portrait of hydrological changes at the 2050 horizon for the province of Québec[END_REF][START_REF] Jennings | Snowfall fraction, cold content, and energy balance changes drive differential response to simulated warming in an alpine and subalpine snowpack[END_REF]). An earlier melt onset implied that net radiation was lower in W20-21, explaining a lower melt rate than in the reference winter (Fig. 5 and Fig. 6). When comparing the subcanopy snow to that found in the gaps, we observed a decrease in melt rate with canopy closure in both years, which is in agreement with previous studies as snowmelt is driven by 𝑆𝑊𝑅 -)3,1( in forest gaps [START_REF] Malle | Shading by trees and fractional snow cover control the subcanopy radiation budget[END_REF]Lawler and Link, 2011;Ellis et al., 2011). However, with less snow accumulating under the canopy due to interception, the snowpack melted-out earlier than in both gaps, in both years. Our results support conclusions from previous studies as canopy interception exerts the main control on snow accumulation patterns whereas distribution of shortwave radiation is the main driver of ablation patterns [START_REF] Lundquist | Lower forest density enhances snow retention in regions with warmer winters: A global framework developed from plot-scale observations and modeling[END_REF]Mazzotti et al., 2023;Dharmadasa et al., 2023). Overall, it shows that the structure of the canopy must be considered in models to precisely predict the effect of climate change on snow accumulation and melt dynamics in boreal forests.

Increased frozen soil under the canopy

A thinner snowpack allows more heat loss from the ground to the atmosphere, and as such, a larger ground heat flux (GHF; Fig. 7). A thinner snowpack in the low-snow winter thus led to a cooler soil than for the reference winter (Fig. 7; Table 5).

This phenomenon was intensified under the canopy, where the snowpack was thinner and the GHF was larger. Our observations show that the heat loss was sufficient to favor soil freezing under the canopy, but not in the forest gaps. These findings are in partial agreement with [START_REF] Stadler | Measurement of frost-induced snowmelt runoff in a forest soil[END_REF], who observed a cooler ground under the canopy, where less snow accumulates. However, the authors also observed ground freezing in a nearby gap. This could be explained by the much lower snow accumulation at their study site compared to what was observed in the MF in both years.

Fig. 8 suggests that the soil has a similar pore space under the canopy and inside gaps, as it saturates at a volumetric water content (VWC) between 0.50 and 0.55 at all monitored sites. However, the seasonal low of VWC is at 0.38, 0.42 and 0.48 in the medium gap, the small gap and under the canopy, respectively, suggesting a higher water retention under the canopy and therefore smaller pores than in the gaps. This implies less potential for water content to increase in response to snowmelt or ROS event. The slower increase in VWC under the canopy at the onset of the 2021 snowmelt (Fig. 8a) also suggests that infiltration was limited, but not completely restricted compared to what was observed in the gaps. Based on previous experiments in the Montmorency Forest, [START_REF] Proulx | Classification of meteorological conditions to assess the potential for concrete frost formation in boreal forest floors[END_REF] concluded that ROS greater than 20 mm followed by a cold spell would cause macropores in the boreal soil to become ice-filled, limiting subsequent infiltration. According to this criterion, the two ROS of December 2020 (42 and 105 mm) should have favored the blockage of soil pores in W20-21. Overall, our results are in accordance with [START_REF] Demand | Influences of macropores on infiltration into seasonally frozen soil[END_REF], who observed a reduced infiltration into frozen soil compared with unfrozen but highly saturated soil during sprinkling experiments over a sandy loam in southern Germany.

Vertical temperature gradients and snow permeability

Our observations show that a thinner snowpack in W20-21 reduced load compaction and therefore limited snow densification, compared to W21-22. A thinner snowpack also led to higher vertical temperature gradients (|𝜕T/𝜕z|) in the low-snow and warm winter (Fig. 9), resulting in more pronounced faceting, lower SSA, and a higher permeability (Ks) than in W21-22 (Fig. 10). In terms of SSA, no significant difference was observed between the two years for the lowermost part of the snowpack, where we observed DH. This is because |𝜕T/𝜕z| in both years was sufficient to favor DH development, and the SSA of DH remained always around 10 ± 2 m 2 kg -1 , regardless of its stage of development [START_REF] Bouchard | Comparison of snowpack structure in gaps and under the canopy in a humid boreal forest[END_REF][START_REF] Domine | Soil moisture, wind speed and depth hoar formation in the Arctic snowpack[END_REF].

In February and March, the |𝜕T/𝜕z|was larger in W20-21 than for the same period the next year. This resulted to a faster development of FC and DH and to a lower SSA in the middle part of the snowpack. Since snow permeability is closely related to snow metamorphism and |𝜕T/𝜕z| (Domine et al., 2013;[START_REF] Domine | Interactions between snow metamorphism and climate: Physical and chemical aspects[END_REF]Calonne et al., 2012), Ks shows a similar pattern to SSA in the snow profiles. When comparing the snowpack under the canopy to the one in both gaps, we observe that temperature gradient metamorphism was enhanced under the canopy, in line with observations from [START_REF] Bouchard | Comparison of snowpack structure in gaps and under the canopy in a humid boreal forest[END_REF].

Interestingly, our field measurements in forest gaps and under a boreal canopy contradict [START_REF] Domine | Interactions between snow metamorphism and climate: Physical and chemical aspects[END_REF]

Reduced spring streamflow

Figure 11 clearly shows that the spring freshet was earlier and reduced in the low-snow and warm winter compared to the reference winter. This suggests that an earlier and slower melt exerted a major influence on streamflow regime. Our results are consistent with the work of [START_REF] Musselman | Slower snowmelt in a warmer world[END_REF], who simulated reduced streamflow runoff with earlier and slower melt.

However, this alone does not explain the large difference in spring runoff between the two years. As shown in Fig. 11, the precipitation for April and May was much lower in 2021 than in the following year (122 mm vs 253 mm). Given a temporal lag of 3 to 5 years between the recharge of the aquifer and groundwater outflow at BEREV-7A [START_REF] Schilling | Quantifying groundwater recharge dynamics and unsaturated zone processes in snow-dominated catchments via on-site dissolved gas analysis[END_REF], snowmelt and precipitation are likely the main contributors to spring streamflow. Overall, although enhanced soil freezing and a higher snowpack permeability in the low-snow year may have contributed to increase spring streamflow in 2021, our results

show that a thin snowpack, combined to dry conditions that spring, significantly reduced the catchment hydrograph as compared to the reference year.

Limitations and shortcomings

In this study, we examined snow accumulation and melt dynamics from highly detailed in situ measurements to establish links between the various processes controlling spring streamflow in discontinuous boreal forest environments. Although this observational study combines snow accumulation and melt dynamics, soil freezing and snow microstructure in one unique dataset representative of the humid boreal forest, the analysis itself remains a case study. Other specific characteristics of the study site such as the slope, the aspect and the surrounding topography, influence the formation and the ablation of the snowpack in forested environments [START_REF] Lundquist | Onset of snowmelt and streamflow in 2004 in the Western United States: How shading may affect spring streamflow timing in a warmer world[END_REF][START_REF] Ellis | Modeling increases in snowmelt yield and desynchronization resulting from forest gap-thinning treatments in a northern mountain headwater basin[END_REF]Mazzotti et al., 2023). To assess the impact of low-snow and warm winter conditions on snowpack dynamics at broad scales, the impact of these factors should be considered.

Some other shortcomings from this study may be relevant for future experimental research projects focusing on the impact of weather conditions on snowpack-related hydrological implications in discontinuous forest at the process scale. One experimental limitation to our study regards the spatial distribution of the soil monitoring. Soil moisture was only monitored at one depth in three nearby locations, which is insufficient to give insights on infiltration patterns at the catchment scale. It underscores the need of monitoring the soil VWC at multiple forest gaps and canopy locations in future studies to better understand soil infiltration dynamics in discontinuous boreal forests under a warming climate. Furthermore, a complete interpretation of soil liquid water content data would benefit from the knowledge of soil granulometry.

Climatic, hydrological and ecological implications

In the low-snow and warm winter, snowmelt started much earlier, but occurred over a longer period at all sites, resulting in an earlier melt-out compared to the reference winter. A shorter-lived snowpack in the boreal forest has climatic impacts, as the net shortwave radiation increases due to a reduction in surface albedo [START_REF] Manninen | Simulation of the effect of snow covered forest floor on the total forest albedo[END_REF] and contributes to a positive feedback loop that enhances global warming [START_REF] Thackeray | Snow albedo feedback: Current knowledge, importance, outstanding issues and future directions[END_REF]. Moreover, a decrease in snow cover extent in the boreal forest may increase the risk of summer hydrological drought due to a lower groundwater recharge (Van Loon et al., 2015). Based on the definition of Van Loon et al. (2015), the conditions in spring 2021 (Fig. 11) have the characteristics of a snowmelt drought.

With increasing warming-induced ROS events, the formation of melt-freeze layers, as we observed at the base of the W20-21 snowpack, and within the snowpack are likely to become more common. These could alter liquid water pathways, favoring water ponding and lateral flow over percolation (Eiriksson et al., 2013;Paquotte and Baraer, 2022;Webb et al., 2018). In addition to altering downward liquid water transport through the snowpack, melt-freeze formations can limit soil-atmosphere gas exchange, promoting hypoxic conditions in subnivean environments [START_REF] Crawford | Oxygen deprivation stress in a changing environment[END_REF]. Melt-freeze formations also limit access to food and movement of subnivean mammals (Johnsen et al., 2017;Poirier et al., 2019), and restrict foraging by large herbivores (Hansen et al., 2011;[START_REF] Schmelzer | Boreal caribou survival in a warming climate, Labrador, Canada 1996-2014[END_REF]. As we observed in December 2020, intense ROS events also trigger winter snowmelt. This promotes low-snow conditions that can intensify soil freezing, as observed under the canopy snowpack in W20-21. Deeper soil freezing in forests decreases microbial activity and soil respiration in winter, thus reducing soil nitrogen recycling (Patel et al., 2018;[START_REF] Yang | Immediate and carry-over effects of increased soil frost on soil respiration and microbial activity in a spruce forest[END_REF]. In turn, this favors carbon accumulation which offsets increased soil respiration in summer (Patel et al., 2018). In summary, it seems clear that changes in snowpack thickness and structure with climate warming have both hydrological and ecological implications in the boreal forest.

Conclusion

Using dedicated field observations, we investigated the effects of a low-snow and warm winter on snow accumulation and melt dynamics, on soil thermal regime and moisture, and on snowpack physical properties under the canopy, and inside two forest gaps in a humid boreal site of eastern Canada. More precisely, we focused on winter 2020-21 (W20-21), which was exceptionally warm and dry, comparing it with a reference winter (W21-22), closer to climate normals. The experimental setup included snowpack and the soil thermal regime monitoring, along with monthly snow pit observations. Our results show that the snowpack was generally half as thick and started to melt earlier in the low-snow and warm winter, compared to the reference winter. This increased soil freezing under the canopy, but not in forest gaps, where the soil remained unfrozen, as for the reference winter. Although the snow surface was warmer in W20-21, the thinner snowpack led to an increase in the vertical temperature gradient, so kinetic growth was enhanced. This resulted in a higher snow permeability during the warm, low-snow year, particularly under the canopy. A thin and early melting snowpack, combined with less precipitation in spring, resulted in a much lower spring runoff in W20-21 than in the reference winter, suggesting that weather conditions were the main driver of spring freshet, rather than soil freezing and snowpack characteristics.

The conditions experienced in the winter 2020-21 at Montmorency Forest, such as warmer air, less snowfall, and a thinner snowpack, were exceptional compared to the past climatology, but may become the new norm in eastern Canada with climate change. Although this work is limited to a two-year comparison within a small catchment, it highlights the many potential effects, all together, of a changing climate on snow hydrology in a discontinuous boreal forest through a unique set of highly detailed process-level observations. These are highly valuable for the snow science community as they will help improve existing modeling tools and develop new ones to address future challenges in snow hydrology.

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

Figure 1 :

 1 Figure 1: Map of the province of Quebec, Canada, with the location of the Montmorency Forest (MF) indicated by a red star (a). Elevation map of the study catchment (BEREV-7A) with the location of the experimental site and the outlet of the catchment (b). Aerial view of the experimental site with the medium gap (yellow), the small gap (purple), and the canopy (green) stations (c) and black and white hemispherical photos of each station location (d-e-f). Picture of the monitoring station under the canopy (g).

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License. classes used for the estimation the reflected shortwave radiation by the snow surface

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License. 2020-21 resulted in the earliest melt-out in the last 40 years at NEIGE site (11 April 2021), with snow disappearing 23 days earlier than the 1982-2022 average (Table

Figure 2 :

 2 Figure 2: Monthly sum of precipitation (a) and means of air temperature (b) from November to May in the winters of 2020-21 and 2021-22, measured at the NEIGE site and compared to the monthly sum of precipitation and means of hourly temperature over the 1982-2022 period. Dashed bars in (a) show snowfall whereas solid bars indicate rainfall. The standard deviation of the 1982-2022 is shown by the black error bars for precipitation in (a) and by the gray area for the monthly temperature in (b). Winter 2020-21 Winter 2021-22 value anomaly rank (out of 40)* value anomaly rank (out of 40)* Ptot JFMA (mm)199 -224 1 428 +5 16 Ptot (Nov-May) 596 -237 2 821 +12 19 Psol (Nov-May) 364 -198 2 596 +13 19 Temp. DJF (°C) -10.7 +2.8 4 -15.2 -1.8 31 Hs,max (cm) 67 -38 5 142 +37 33 Melt-out day (DOY) 101 -23 1 134 +11 33 * Temp. DJF is ranked in descending order, whereas all the other variables are ranked in ascending order.Table 4: Total precipitation from January to April (JFMA) and from November to May, solid precipitation from November to May, mean December-January-February (DJF) temperature, maximum recorded snow height (Hs,max), and the melt-out day at the NEIGE site for the 2020-21 and 2021-22 winters. The anomaly relative to the 1982-2022 period is also shown, along with the rank of both winters for each metric relative to the 40 winters of the climatology. We used a threshold at 1°C to define the precipitation phase over the analysis period. Winters 1999 to 2003 were excluded of the analysis because of too many missing data.

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

Figure 3 :

 3 Figure 3: Snowpack thermal regime and cold content in the medium gap (a-b), in the small gap (c-d) and under the canopy (e-f) for winters 2020-21 (left) and 2021-22 (right). Dashed lines show the onset of snowmelt.265

Figure 4 :

 4 Figure 4: Snow surface and air temperature measured in the medium gap (yellow), in the small gap (purple) and under the canopy (green) during winters 2020-21 (a-c) and 2021-22 (b-d). Since the snow surface temperature in the small gap was entirely gap-filled in winter 21-22, it is not presented in (b).

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

Figure 5 :

 5 Figure 5: Cumulative net total (Rnet,bc), shortwave (SWRnet,bc) and longwave (LWRnet,bc) radiation below the canopy during snowmelt 280

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

Figure 6 :

 6 Figure 6: SWE evolution during snowmelt of winters 2020-21 (red) and 2021-22 (blue) in the medium gap (a), in the small gap (b) and under the canopy (c). The plots start at the beginning of the snowmelt period on each year, which differs by an average of 23 days between seasons.

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

Figure 7 :

 7 Figure 7: Ground heat flux (GHF) and temperature at 20, 10, 5 cm below ground level and at the soil-snow interface during the snow cover period of winters 2020-21 (left) and 2021-22 (right) in the medium gap (a-b), in the small gap (c-d) and under the canopy (e-f).305

Figure 8 :

 8 Figure 8: Soil liquid water content (VWC) at 15 cm below the surface in the medium gap (yellow), in the small gap (purple) and under the canopy (green) in winters 2020-21 (a) and 2021-22 (b).

Figure 9 :

 9 Figure 9: 7 day rolling mean of the magnitude of the vertical temperature gradient (|𝝏T/𝝏z|) in winters 2020-21 (red) and 2021-22 (blue) in the medium gap (a), the small gap (b) and under the canopy (c). The grey band on each frame shows the transition zone from equilibrium to kinetic growth.

Figure 10

 10 Figure 10 presents the comparison of snow stratigraphy, density, SSA and permeability between both years for gaps and under the canopy as obtained from snow pit observations between 9 and 11 March both years. These snow pit dates are convenient to present as they are similar in both years and correspond to a snowpack under dry snow conditions and well into the snow accumulation period. Observations from the other snow pits are presented in Supplementary Material (Fig S8 -S11).

Figure 10 :

 10 Figure 10: Stratigraphy, snow density, SSA, and snow permeability profiles as measured in the medium gap (a-b-c-d), in the small gap (e-f-g-h) and under the canopy (i-j-k-l) on 9 and 10 March 2021 (red) and 10 and 11 March 2022 (blue).

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

Figure 12 :

 12 Figure 12: Summary of the results. Upward arrows correspond to an increase and downward arrows to a decrease in the low-snow and warm winter with respect to the reference winter. The clock with counterclockwise arrows means that the process happens earlier. The yellow, purple, and green arrows indicate the effects in the medium gap, the small gap and under the canopy, respectively. Their relative size shows the importance of the process between the gaps and the subcanopy locations. Large black arrows are applied all three locations. Small black arrows show the causal link between the observations processes. Gray boxes refer to processes treated in this study (1. snowmelt dynamics; 2. soil thermal regime; 3. snow metamorphism).

  -e and Fig. S3.1a. Based on these weather and snowpack https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License. conditions, we consider the winter 2020-21 to be representative of what could be expected for eastern Canada with global warming.

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

  https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

Table 1 : Gap dimensions and sky-view fraction of each station

 1

Table 2 : Snow pit measurement dates

 2

		Medium gap	Small gap	Canopy
		8 Dec. 20	8 Dec. 20	8 Dec. 20
	Winter 2020-21	27 Jan. 21	27 Jan. 21	26 Jan. 21
	(low-snow)	9 Mar. 21	10 Mar. 21	10 Mar. 21
		6 Apr. 21	6 Apr. 21	6 Apr. 21
		13 Jan. 22	12 Jan. 22	12 Jan. 22
	Winter 2021-22	15 Feb. 22	15 Feb. 22	15 Feb. 22
	(reference)	10 Mar. 22	11 Mar. 22	10 Mar. 22
		28 Apr. 22	29 Apr. 22	28 Apr. 22

https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

  who, based on general considerations, suggested a decrease in snow permeability with global warming mainly due to increased snow surface temperature. Our results suggest that if warmer winters lead to a large decrease in snowpack thickness, this could override the increase in surface temperature and lead to higher |𝜕T/𝜕z| and permeability.

https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

https://doi.org/10.5194/hess-2023-191 Preprint. Discussion started: 31 August 2023 c Author(s) 2023. CC BY 4.0 License.

Figure 11: Air temperature (a-b), daily difference of SWE, total precipitation (c-d) and streamflow discharge (e-f) for winter 20-21

Acknowledgements. The authors thank the staff of the Montmorency Forest for helping us logistically with field visits. We also thank Charles Villeneuve and Kino Leroux for preparing snowmobiles trails before our visits. The authors also thank Éric Boucher, Christian Juneau and Antoine Thiboult for helping in the deployment and the maintain of the monitoring stations.

We thank Pierre-Erik Isabelle for providing radiation dataset. We finally thank all the people that accompanied Benjamin Bouchard and Etienne Tremblay on the field, especially the members of PÉGEAUX and the graduate students, post-doc, and research associates from the Hydrometeorology Lab at Laval University. The work of BB was founded by the Natural Sciences and Engineering Research Council (NSERC) and Sentinel North program. The authors acknowledge the financial contribution of Environment and Climate Change Canada through the Grants & Contributions program (projects #GCXE20M016 and #GCXE22M013) and the Cold-region climate projection for hydrological applications (EVAP-2; project #ALLRP 549108 -19).

Data availability. Data from snow pit measurements and from the monitoring stations in the medium and small gaps, and under the canopy are freely available at https://doi.org/10.5281/zenodo.8213204.

Author contributions. BB, DFN and FD designed the study. BB and ET collected and treated field data. BB conducted the analysis of the results with inputs from DFN, FD, FA, and TJ. BB wrote the manuscript, with feedbacks from all authors.

Competing interests. The authors declare that they have no conflict of interest.


Smart Nanomaterials

"Smart" materials-materials that respond to a stimulus or their environment to produce a dynamic and reversible change in critical properties-have enabled progress in many areas, including display technologies, drug delivery, and selfhealing materials for coating applications, among others.

ACS Nano 2008Nano , 2, 6, 1101Nano -1107 Magnetic and Magneto-Optical Nanoparticles J. Chem. Theory Comput. 2022, 18, 5, 3122-3135 Plasmon induced magneto-optical enhancement in metallic Ag/FeCo core/shell nanoparticles synthesized by colloidal chemistry Nanoscale, 2018,10, 18672-18679

Zeeman Alignment and Dipolar Assembly of Magnetic Nanoparticles

The combination of magnetic and plasmonic materials and their nanostructuring have revealed a prominent pathway to develop novel photonic materials for the active control of the light polarization using a magnetic field. routine that corrects n, k for small particles (mean free path) Haiss et al. Analytical Chemistry, Vol. 79, No. 11 (2007) ok for size NP >5nm

  propagation, electronic excitation, light scattering, reflection and absorption • Electron-phonon relaxation, heat diffusion, pressure wave propagation, rarefaction and tensile stress, phase transitions, mechanical fragmentation • Materal removal (ablation), plasma plume expansion, reactions, clustering, etc… • Often multi-pulse irradiation is used with various scanning and focusing conditions Physics of Laser Interactions Laser-induced nanoparticle formation • Green synthesis methods, biocompatibility • Control over NP properties • Scaling-up possibilities • Possibility to create stable ultrastructures, nanoalloys, hybrids, etc. • Additional laser-particle interactions involve many physical and chemical processes In gas or in vacuum with different geometry In liquids or in the precence of a Process can be at T-> Tc Difficult to model by classical methods pico T. E. Itina, L. V. Zhigilei, B. J. Garrison, J. Phys. Chem.B , 106, 303-310 (2002) Short (ps) and ultra-short (fs) laser pulses ps, F ~ 2 Fth, 4010 nm, t =1 ns Collaboration with L.V. Zhigilei et al. (Virginia State University, USA) Practically all the molten material is ablated Labl ~ L0 ln(Fabs/F0), E>Eb . Itina, L. V. Zhigilei, B. J. Garrison The Journal of Physical Chemistry B 106 (2), 303-310 (2002) Laser Interactions with Nanoparticles Laser pulse Propagation, absorption, local field enhancement Local ionization and heating of the matter/liquid near NP plasma-mediated nanobubbles NP heating, heat transfer to the ambient medium thermo-mediated nanobubbles NP heating, melting, reshaping, sintering, evaporation Electron ejection, injection in the ambient NP fragmentation Lachaine et al.Photonics (2013) Hashimoto et al. J. PhotChem&PhotBiology (2012) Strasser et al, . J. Phys.Chem. C 118, 25748 (the medium. At 400 nm, passage from liquid water to water vapor goes with a drop in absorption. Absorption at 400 nm in liquid water raises with the particle radius (not linear) , absorption cross section σ ext , σ scat , σ abs A generalized multiparticle Mie-solution: further experimental verification Yu-lin Xu et al. FDTD, EM FEM, DDSCAT, etc..

+

  Bohren-Huffman Mie scattering subroutine (callbhmie.f) Generalized Multiparticle Mie-solution (gmm01f.f) based on

Bimetal CuFe Nanoparticles: Synthesis and Electro-catalytic Activity Electro-Catalytic Activity of CuFe Nanoparticles

Zaneta

Swiatkowska-Warkocka, Appl. Sci. 2021, 11(5), 1978

Rojas-Nunez, et al.Toward Controlled Morphology of FeCu Nanoparticles: Cu Concentration and Size Effects. J.

Phys. Chem. C 2018, 122, 8528-8534.

Ultra-short Laser-induced NP sintering, Janus NPs Sintering of a larger metallic Janus NPs with different sizes upon rather weak laser heating

nm radius

tatiana.itina@univ-st-etienne.fr

Thank you very much!

Just published:

"Alternative local melting-solidification of suspended nanoparticles for heterostructure formation enabled by pulsed laser irradiation" Advanced Functional Materials 2023.


Lorenzo Cabona

Edoardo Redaelli

Frederic Zamkotsian

Giorgio Pariani

Patrick Lanzoni

Andrea Bianco

Silicon based photo-controlled deformable mirror

Keywords: Deformable mirror, Photo-controlled deformable mirror, Wavefront shaping, Adaptive Optics

come

Introduction

Deformable mirrors are optical devices that are finding an increasing interest in a wide range of applications including wavefront correction [START_REF] Tyson | Adaptive optics engineering handbook[END_REF], temporal and spatial beam-shaping [START_REF] Dalimier | Comparative analysis of deformable mirrors for ocular adaptive optics[END_REF], and active and adaptive optics [START_REF] Brida | Sub-twocycle light pulses at 1.6 $µ$m from an optical parametric amplifier[END_REF]. In the case of astronomical telescopes, they allow for correcting aberrated wavefronts (by the atmosphere or the inhomogeneous medium) and reach a better resolution of the optical system [START_REF] Rigaut | Multiconjugate adaptive optics for astronomy[END_REF]. This is especially true for large size telescopes that have a very high theoretical resolution. This is possible thanks to the capability of such devices to modify their optical surface by suitable actuating systems [START_REF] Madec | Overview of deformable mirror technologies for adaptive optics and astronomy[END_REF]. Conventional technologies for the decimetres DMs are mainly based on two actuation processes: electromagnetic for the telescope adaptive secondary mirrors, piezoelectric for the DMs in the AO modules; the first ones are built around voice coils and thin shell mirrors [START_REF] Madec | Overview of deformable mirror technologies for adaptive optics and astronomy[END_REF] [START_REF] Biasi | Contactless thin adaptive mirror technology: past, present, and future[END_REF], the second ones use stack actuators [START_REF] Wlodarczyk | Scalable stacked array piezoelectric deformable mirror for astronomy and laser processing applications[END_REF], or annular PZT (Lead Zirconate Titanate) layers for unimorph bimorph DMs [START_REF] Zhu | Development of a unimorph deformable mirror with water cooling[END_REF]. The development of new technologies based on microopto-electro-mechanical systems (MOEMS) is promising for future centimetre-size DMs. Their major advantages are compactness, scalability, and specific task customization using elementary building blocks [9]. The MOEMS technology permits the development of a complete generation of new mirrors. Four main Micro-Deformable Mirrors (MDM) architectures are under study in different laboratories and companies, and they are based on single-crystalline and poly-crystalline silicon materials for defining the structural layers, i.e. mirrors, restoring beams, and hinges. First, the bulk micro-machined continuous-membrane deformable mirror, studied by Delft University and OKO company, is a combination of bulk silicon micromachining with standard electronics technology [START_REF] Vdovin | Deformable mirror with thermal actuators[END_REF]. This mirror is composed by a thin flexible conducting membrane, coated with a reflective material, and stretched over an electrostatic electrode structure. This mirror shows a very good mirror quality, but the mean deformed surface is concave, and the number of actuators cannot be scaled up to hundreds. Second, the segmented, microelectro-mechanical deformable mirror manufactured by Iris AO [START_REF] Helmbrecht | Long-term stability and temperature variability of Iris AO segmented MEMS deformable mirrors[END_REF] consists of a set of segmented piston-tip-tilt moving surfaces, fabricated in a dense array. Third, the surface micro-machined continuous-membrane deformable mirror made by Boston Micromachines Corporation (BMC) is based on a single compliant optical membrane supported by multiple attachments to an underlying array of surface-normal electrostatic actuators [START_REF] Cornelissen | MEMS deformable mirror actuators with enhanced reliability[END_REF]. The efficiency of this device has been demonstrated recently in several AO system for astronomy, including the GPI instrument on Gemini telescope [START_REF] Macintosh | The Gemini planet imager[END_REF]. The third concept is certainly the most promising architecture, but it shows limited strokes for large driving voltages, and mirror surface quality may need further improvement for Extreme AO. Another interesting technology is the one employed in the ALPAO DMs [START_REF] Vidal | Tests and characterisations of the ALPAO 64 x 64 deformable mirror, the MICADO-MAORY SCAO AIT facility[END_REF], that are based on continuous reflective surface deformed by magnetic actuators, in this case the maximum stroke for tip/tilt reaches 90 .

In all the above-mentioned technologies the density of actuators cannot be increased independently of the technology. Besides these technologies a deformable mirror can be realized using a large and continuous metalized membrane and a monolithic nonpixelated photoconductive substrate with the possibility to avoid the presence of real, physical actuators. Such devices are called Photo Controlled Deformable Mirrors (PCDMs).

The PCDM is an electrostatic deformable mirror, composed by a thin reflective membrane suspended above and electrically coupled to a photoconductive slab [15] [16]. When the photoconductive slab is illuminated with a light pattern, a change in the electrostatic pressure is produced in the illuminated areas of the photoconductor, which locally modifies the shape of the membrane (see Figure 1). Appropriately managing the light pattern illuminating the photoconductor, it is possible to define the size, shape and density of actuators, without the spatial segmentation into discrete elements inevitably present using the other deformable mirrors technologies [START_REF] Bortolozzo | Continuous photocontrolled deformable membrane mirror[END_REF]. The particular design of the PCDM takes advantage of well-developed light projection techniques, allowing to decouple the driving electronics from the mirror actuation. This may be a huge advantage in comparison to usual AO mirrors that have a bulky and complex electronics. This turns into a potential reduction of cost, weight and space occupied by mechanical actuators [START_REF] Bonora | High resolution wavefront correction with photocontrolled deformable mirror[END_REF]. These features make the PCDM technology very attractive in those situations where the weight and volume of the deformable mirror must be reduced [18] [19].

PCDMs have been manufactured using crystalline inorganic photoconductors such as Gallium Arsenide (GaAs), Indium phosphide (InP) [START_REF] Haji-Saeed | Photoconductive optically driven deformable membrane for spatial light modulator applications utilizing GaAs and InP substrates[END_REF], Bismuth-Silicon-Oxyde (BSO) [START_REF] Bonora | High resolution wavefront correction with photocontrolled deformable mirror[END_REF][15], Zinc Selenide (ZnSe) [START_REF] Quintavalla | Zinc selenide-based large aperture photo-controlled deformable mirror[END_REF] and, more recently, using an organic photoconductor consisting in a mixture of P3HT/PCBM doped with an arylamine (m-MTDATA) dispersed in a polystyrene matrix [START_REF] Quintavalla | Fully Organic Photocontrolled Deformable Mirror[END_REF]. A model which identifies the key parameters that affect their performances is reported in the literature [START_REF] Bortolozzo | Continuous photocontrolled deformable membrane mirror[END_REF]. To the best of our knowledge, silicon has not been used so far for building a PCDM. This material is attractive for the high charge mobility; moreover, it is easy to find at a relatively low price and in large formats.

In this paper, we show the manufacturing of a 1-inch silicon based PCDM and its characterization in terms of deformation as function of the driving parameters. In addition, the response time is evaluated.

Material and methods

Photoconductor

The photoconductor is a round 1-inch undoped silicon substrate manufactured by UniversityWafer, Inc. The main properties of the photoconductor slab are reported in Table 1. On one side of the Silicon slab, a conductive ITO coating has been deposited by sputtering followed by an annealing at 400°C for one hour. The thickness of the ITO layer is in the 100 nm range and the final resistance is 100 Ω/sq. The transmittance of the bare and coated silicon substrates is reported in Figure 2. As expected, the transmittance increases when the ITO coating is applied thanks to the reduction of the reflection losses.

Reflective membrane

The membrane used in our PCDM were purchased from National Photocolor. They are a 5 µm thick nitrocellulose layer with a reflective Al coating. The membrane is mounted on a rigid aluminium ring with a 25.4 mm internal diameter and 35 mm outside diameter. The membranes were characterized in terms of wavefront without voltage applied and tension. These tasks were carried out using an aluminium dummy photoconductor (uniform load) and a Zygo GPI interferometer. The wavefront of the bare membrane is reported in Figure 3. The wavefront is astigmatic probably due to the gluing of the membrane to the supporting ring. The astigmatism of the membrane introduces a Peak-to-Valley (PV) deformation of 940 nm. To evaluate the membrane tension, a DC voltage between 0 to 180V were applied in ten steps to a dummy aluminium disk electrically coupled with the membrane. Its maximum displacement was measured at the centre as a function of the applied electric field and the data fit with the thin membrane equation: [START_REF] Tyson | Adaptive optics engineering handbook[END_REF] where ε is the dielectric constant, φ the membrane diameter, V the supply voltage, d the distance between the dummy photoconductor and the membrane and M the maximum static deformation. From the tension values, the first resonance frequency in vacuum was also calculated make use of the formula:

. ( 2
)
where σ is the mass/surface ratio (for nitrocellulose σ is 1230 / ). The results for the employed membrane are summarized in Table 2.

Mechanical design, Interfaces and assembling

The design of the PCDM device needs to guarantee some mechanical specifications, in particular:

• Guarantee the stability in position of the components during the installation procedure; • Keep the charged surfaces (photoconductor and reflective membrane) well isolated between each other; • Guarantee the possibility to align and compensate the manufacturing tolerance for both the photoconductor and the membrane; • Keep the mechanics centered and able to be aligned with respect to the external optical device.

For this reason, the adopted strategy was to use calibrated shims between the photoconductor and the membrane supports. These calibrated shims are able to guarantee a tip/tilt resolution of 0.1 mrad. A key part of the PCDM device was the holder for the membrane and the photoconductor that must be made of a non-conductive material. We opted for acrylonitrile butadiene styrene (ABS). The mounting allows for shining the photoconductor from the back. More details on the mechanical aspects are reported in a previous publication [START_REF] Redaelli | Photocontrolled deformable mirror based on silicon photoconductor[END_REF]. Figure 4 shows the 3D printed components of the PCDM and the exploded view drawing of the assembled PCDM.

Characterization setup

The characterization of the PCDM were performed measuring the wavefront of a collimated He-Ne laser beam reflected by the membrane. The wavefront was recorded by means of the PhotonLoop Shack-Hartman wavefront sensor [SHWS]. It is equipped with an iDS UI-3060 CMOS camera controlled by the PhotonProbe software [START_REF] Mocci | A multiplatform CPUbased architecture for cost-effective adaptive optics systems[END_REF]. The main features of the Shack-Hartman device are reported in Table 3. The scheme of the optical set up is reported in Figure 5 and described in [START_REF] Zamkotsian | Photocontrolled deformable mirror for wavefront shaping[END_REF]. The PCDM was fed with an AC voltage obtained by the combination of a function generator (Voltcraft FG8210) and a voltage amplifier (Falco Systems WMA-100). For the characterization measurements the range of the applied voltage was from -180 up to 180 V and the frequency from 3 kHz to 30 kHz. The illumination of the devices was performed on the back side (surface with the conductive ITO coating) of the silicon photoconductor with a NIR LED or lasers. The LED (Roithner LaserTechnik smb1n-940d-02) emits light peaked at 940 nm and provides an irradiance in the range of 2.8 -8.45 / on the ITO coating. The two lasers have a wavelength of 532 nm and 1064 nm, the first with a power of 110mW and the latter of 11mW. The spot diameter is 2.5 mm. The switching on and off of the LED and the duty cycle was achieved by means of a microcontroller board (Arduino Mega 2560). The time response of the PCDM was evaluated using a fast photodiode as a detector. The photodiode is a Texas Instruments OPT101P, with a bandwidth of 14 kHz and a rise and fall time (10%-90%) of 28 µs. To perform this measurement, the 940 nm LED was used to illuminate the PCDM. The LED showed a response time < 100 µs. When the light is illuminating the photoconductor of the PCDM, the curvature of the membrane changes; in order to discriminate between the ON and OFF states of the device, we locate a filtering hole in an intermediate focal plane. By turning on the illumination a drop of the analog output of the signal generated by the photodiode was detected. The temporal response of the PCDM could then be recorded.

Results and discussion

The PCDM was firstly characterized in terms of maximum deformation as function of the driving parameters. The photoconductor was illuminated on the back side with a uniform light at 940 nm as described in the previous section. Figure 6 shows the interference fringes for the membrane measured with the SHWS under various conditions of applied voltage and illumination. Figure 6a shows the interference fringes when a voltage of 120 V at 20 kHz is applied, in dark conditions and without absolute reference, Figure 6b shows the interference fringes in the same conditions but with a light intensity of 5.63 / , Figure 6c shows the difference between cases b) and a). In order to characterize the PCDM, the SHWS was set to zero in dark condition with the AC voltage applied; therefore, once the light was turned on, the measured deformation was only due to the illumination. The run consisted in collecting the wavefront sampled at 10 Hz in 20 on-off cycles (1 second with the light on and 1 second with light off). After we have acquired the wavefront, the mirror surface is reconstructed using the first 37 orders of the Zernike polynomials calculated and recorded by the PhotonProbe software of our SHWS. Then, we were able to retrieve the maximum deformation of the membrane, i.e., the difference in Z from the edge to the center of the membrane, with respect to time (see Figure 7a). Finally, for each on-off cycle was computed the mean value and the standard deviation of the maximum membrane deformation of each on-off cycle.

In the bottom of Figure 7, the trends of this parameter with the applied voltage, its frequency and light irradiance are reported.

Figure 7: An example of the measurements performed to retrieve the plots with 20 cycles 1s ON + 1s OFF (a). Root mean square (rms) of the maximum deformation of the membrane and their fit (red lines) as function of the applied voltage (AC, peak) (b), voltage frequency (sinusoidal) (c), and irradiance at 940 nm (d). The bars represent the standard deviations of the rms with respect to each set of 20 on-off cycles. The mechanical deformation, in the parameter's interval considered, shows a square dependence with the applied voltage, a linear dependence with the voltage frequency and a negative exponential dependence with the Optical Power.

We notice an increase of the deformation with the applied voltage and this is consistent with the working mechanism of the PCDM and in agreement with the results reported on similar devices [START_REF] Quintavalla | Zinc selenide-based large aperture photo-controlled deformable mirror[END_REF]. The deformation also increases with the AC frequency and this can be explained considering that the resistance of the silicon is not very high, but it becomes less and less important when increasing the frequency with respect to the capacitor impedance. It is important to stress the fact that the frequencies are well above the resonance frequency of the membranes. We have an increase of the deformation with the LED irradiance and the trend is towards a plateau, probably due to the saturation effect of the photovoltage responsivity, and photocurrent responsivity [START_REF] Paschotta | Encyclopedia of Laser Physics and Technology -responsivity, photodetectors, photodiodes, sensitivity[END_REF]. Again, this is consistent with the charge generation and transport and it has been found in other PCDMs [21][22].

Concerning the values, we obtain a maximum mechanical deformation in the center of the membrane up to 1.4 µm, this value is below the deformation obtained in similar membrane DMs, such as OKO which reach deformations of the order of 10 µm [10], but comparable with the surface micro-machined continuous-membrane DMs, such as Boston Micromachines. We have to consider that, in PCDM, the applied voltage induces a constant deformation, which is not present in the common electrostatic membrane mirrors.

Another important task was the determination of the temporal stability of the induced deformation. This was evaluated performing a sequence of 30 seconds light ON and 30 second light OFF and calculating the stability of the maximum deformation. In Figure 8, we report the results for two datasets. The results are comparable in terms of achieved deformation. Concerning the stability, the dataset A shows a rms value of 0.022 µm and a PV of 0.084 µm in the ON state and a rms value of 0.012 µm and a PV of 0.039 µm in the OFF state, while for the dataset B shows a rms value of 0.010 µm and an PV of 0.050 µm in the ON state and a rms value of 0.007 µm and a PV of 0.041 µm in the OFF state. The first dataset is nosier in the ON state than the second, which shows similar values for the ON and OFF states. This is linked to a relative instability of the the location of the maximum deformation, probably due to the environmental factors (e.g., vibrations, air turbulence and light intensity variation). This time scale is sufficient for applications in adaptive optics where the refresh rate is of the order of 0.1 -1 kHz.

For testing the long-term stability of the device, a set of 220 1s on -1s off repetitions was performed. The results are reported in Figure 9 in terms of maximum deformation. We see the very high stability of the undeformed position (light OFF) consistent with the previous results and a noisier value for the ON state with a slight drift towards smaller deformations (the maximum mechanical deformation was decreased less than 8%). Indeed, the average value for the OFF and ON states are 0.0488 ± 0.0212 and 1.0297 ± 0.0379 respectively. This outcome is not so critical in the case of closed loop adaptive optics system, where the deformation is driven at a high rate. After the tests with the uniform illumination that showed the general performances of the PCDM, lasers at 532 nm (110mW) and 1064 nm (11mW) were projected on the back side of the device in order to probe the deformation with a confined virtual actuator. The spot size was 2.5 mm in diameter and it was scanned along the wafer radius. The membrane deformation was evaluated following the established protocol. In Figure 10, the deformation profile is reported for the different spot positions. We notice how the deformation is well localized in the illuminated area, but it is spread over almost half of the membrane area, illustrating a very large influence function. Considering the Full Width at Half Maximum (FWHM), it is about 12 mm, which is half of the PCDM diameter. This could be due to the spread of the charges in the horizontal direction after the photogeneration along the photoconductor substrate, due to the relative high mobility in this direction promoting the charge diffusion and the fact that the sinusoidal vertical electric field is not very high. Indeed, at 20 kHz the voltage drop across the photoconductor is about 12% of the applied voltage, meaning that for an AC 120 V, only 10 V drops on the photoconductor. We can than locate about 10 lasers spots on the surface, i.e., ten actuators, with a potential correction of focus and other low orders aberrations. Some strategies could be possible to improve the confinement and to narrow the influence functions with respect to the illumination pattern dimensions. The employed lasers at 532 nm and 1064 nm to generate a confined virtual actuator have a higher power than the extended LEDs, but as shown in Figure 7d, we reached the plateau value of the deformation according to the exponential trend with these powers. Therefore, we can induce a similar deformation with much smaller intensities of the order of some mW/cm 2 . The last test concerned the time response of the device and it has been evaluated following the procedure reported in the Materials and Methods section. In Figure 11, the signal of the photodiode as a function of time is reported. The intensity drop-rise measured by the photodiode is sudden after the light OFF-ON trigger on the PCDM. The temporal response was characterized by illuminating the PCDM locally on an area of 1 , simulating the action of a single conventional actuator and measuring the membrane shape variation when the light stimulus is applied (rise time) and when light is removed (fall time). The time response of the device was quantified fitting each single cycle with a simple negative exponential function both for the rise and fall and then averaging the time constant τ parameter. The average τ is equal to 0.0052 ± 0.0003 s for the fall and 0.0065 ± 0.0006 for the rise. Considering the 10% to 90% deformation (it is equal to 2.2 τ), we obtain 0.0115 ± 0.0005 s and 0.0142 ± 0.0012 s respectively. We notice a slighter slower process when the light is turned OFF and this is consistent with other results reported in the literature [START_REF] Quintavalla | Zinc selenide-based large aperture photo-controlled deformable mirror[END_REF] and can be ascribed to the force that is opposite to the applied pressure. We can conclude that the average frequency for the full deformation range is of the order of 100 Hz.

Conclusions

We have proposed and fabricated a photo-controlled deformable mirror (PCDM) based on an intrinsic 1-inch silicon wafer. A dedicated optomechanical design was performed in order to control the key parameters that affect the PCDM response. The device worked in the AC regime and the deformation was induced by illuminating it with both a NIR LED (940 nm) and lasers (532 nm, 1064 nm). The maximum deformation induced by a uniform illumination was measured by a SHWF and it was a function of the driving voltage, its frequency and the light irradiance, which was the most effective parameter. The largest value was about 2 µm, which is comparable with other deformable mirrors and it allows for a fine control of the mirror focus. By the laser measurements, we identified a low confinement of the deformation, which turned into large influence function that makes the device suitable only for low orders correction of the wavefront. The response time of the Si based PCDM was of the order of 10 ms for the 10%-90% deformation, which turns into a working frequency of 100 Hz.

Figure 1 :

 1 Figure 1: scheme of the functioning of a PCDM. (a) The deformation of the membrane due to the applied tension in dark condition; (b) deformation of the membrane when the photoconductor is illuminated by the light beam, in this last situation the deformation due to the photogenerated electron-hole pairs is superimposed on the previous deformation.

Figure 2 :

 2 Figure 2: transmission spectra of the silicon photoconductor uncoated (blue curve) and coated (orange curve).

Figure 3 :

 3 Figure 3: the interference fringe map of the membrane, see the interference fringes for the mounted membrane without any applied electric field, without voltage being applied and in dark conditions.

Figure 4 :

 4 Figure 4: images of the back side (a) and lateral side (b) of the 3D printed components, image of the assembled PCDM in its optomechanics (c), and exploded view drawing of the assembled PCDM design (d), where the listed components are: 1. Screw with membrane side Al ring, 2. Main Al cylinder support, 3. Photoconductor and ABS Photoconductor support, 4 Separating shims, 5. Membrane ABS holder, 6. Membrane, 7. Plastic cover body, 8. Screw with photoconductor side Al ring, and 9. Mechanical adapter holder.

Figure 5

 5 Figure 5 Photo a) and scheme b) of the measurement setup with the main components: 1) amplifier, 2) function generator, 3) DC power supply, 4) laser, 5) PC , 6) Shack-Hartmann, 7) PCDM and PCDM support, and 8) Light surce and collimation lens.

Figure 6 :

 6 Figure 6: interference fringes of the membrane #2 with an applied voltage of 120 V at 20KHz: in dark condition without absolute reference(a), with a light intensity of 5.63 mW/cm^2 without absolute reference (b), with a light intensity of 5.63 mW/cm^2 after had acquired as absolute reference the condition a (c).

Figure 8 :

 8 Figure 8: maximum deformation as function of time with 30 seconds ON-OFF cycles for the two datasets (the blue line shows the dataset A, while the orange the dataset B). The PCDM is driven at 120 V, 20 kHz and light intensity of 5.63 / 2 .

Figure 9 :

 9 Figure 9: maximum deformation for a set of 220 repetition on-off (2 seconds total period). The PCDM is driven at 120 V, 20 kHz. Light intensity of 5.63 / .

Figure 10 :

 10 Figure 10: on top: plot profiles of the deformed membrane for three different illuminations with laser spots (532 nm, 1064 nm). On the bottom: Scheme of the PCDM reporting the position where the laser spots illuminated the photoconductor (the blue beam shows the position 1, while the one orange the position 2 and the yellow the position 3).

Figure 11 :

 11 Figure 11: photodiode signal with the DMD always on and the led switch on-off by the Arduino bord (a), fits of the membrane maximum mechanical deformation during the lighting on (red) and lighting off (yellow), The PCDM is driven at 120 V,10 Hz. Light intensity of 5.63 / at 940 nm.

Table 1 :

 1 properties of the Silicon photoconductor slab employed in the PCDM

	Property	value
	Diameter (mm)	25.4
	Thickness (mm)	3
	Capacitance (pF)	13
	Structure	Single Crystal
	Resistivity (Ohm cm)	>10 4

Table 2 :

 2 properties of the membrane employed in the PCDM

	Property	Value
	Diameter (mm)	25.4
	Tension [N/m]	77
	Resonance frequency [Hz]* σ [ / ]	3370 1230

* It is calculated considering the membrane in vacuum

Table 3 :

 3 Main features of the Shack-Hartman wavefront sensor.

	Property	value
	Pixel sensor	1535x1216
	Pixel size (µm)	5.86
	Lenslet pitch (µm)	150
	Focal length (mm)	5.2
	Optical area (mm × mm)	9x7.1
	Max number of centroids	48x60
	Max framerate	100Hz with 1200 centroids

Acknowledgements

We are grateful to Dario Natali for the fruitful discussion on the electric behaviour of the device. This work was partially supported by OPTICON project (funded by the European Union's Horizon 2020 research and innovation programme under grant agreement No 730890)