Estimating Snowpack Density from Near-Infrared Spectral Reﬂectance Using a Hybrid Model

: Improving the estimation of snow density is a key task in current snow research. Char-acterization of the variability of density in time and space is essential for the estimation of water equivalent, hydroelectric power production, assessment of natural hazards (avalanches, ﬂoods, etc.). Hyperspectral imaging is proving to be a promising and reliable tool for monitoring and estimating this physical property. Indeed, the spectral reﬂectance of snow is partly controlled by changes in its physical properties, particularly in the near-infrared (NIR) part of the spectrum. For this purpose, several models have been designed to estimate snow density from spectral information. However, none has yet achieved signiﬁcant performance. One of the major difﬁculties is that the relationship between snow density and spectral reﬂectance is non-bijective (surjective). Indeed, several reﬂectance amplitudes can be associated with the same density and vice versa, so the correlation between density and spectral reﬂectance can be very poor. To resolve this issue, a hybrid snow density estimation model based on spectral data is proposed in this work. The principle behind this model is to classify the snow density prior to its estimation by means of a speciﬁc estimator corresponding to a predeter-mined snow density class. These additional steps eliminate the surjective relation by converting it into three bijective relations between density and spectral reﬂectance. The calibration step showed that the densities included within the three classes are sensitive to different spectral regions, with R 2 > 0.80. The results of the cross-validation for the speciﬁc estimators were also satisfactory with R 2 > 0.78 and RMSE < 36.36 kg m − 3 . The overall performance of the hybrid model (HM), when tested with independent data, demonstrated the effectiveness of using proximal NIR hyperspectral imagery to estimate snow density (R 2 = NASH = 0.93).


Introduction
The spatiotemporal evolution of seasonal snowpacks is an important indicator of climate [1]. The measurement, monitoring, and management of this water resource are of great interest to governments and the scientific community. Snow cover is the set of snow layers that accumulate on the ground throughout the winter [2]. Each of these layers has a given density. Quantifying the variability of density in time and space is essential for estimating the water equivalent [3,4], hydroelectric power production [5,6], and assessing natural hazards (avalanches, floods, etc.) [7,8]. Indeed, this variable varies with changes in other physical properties such as grain size, grain shape, and liquid water content during the metamorphic transformation of the snowpack [9]. The density of newly deposited snow is expected to have the lowest values and to increase and reach its highest values during the maturation phase. According to Pomeroy et al. [10], the typical seasonal density of the snowpack ranges between 80 kg m −3 and 600 kg m −3 .
For a country like Canada, which covers a very large area with a vast expanse of snow cover, regular monitoring of snow density is important [11,12]. The density is measured using a variety of methods and technologies. These include manual measurements by taking core samples from the snowpack (such as 'federal' snow tubes, e.g., ESC-30) [13], or the installation of devices that lie flat on the ground and weigh the snow as it accumulates on top (such as snow pillows) [14,15]. However, each of these methods has several drawbacks [15,16]. Snow core measurements are labor intensive, time-consuming, not feasible for 24-h data collection, and subject to human error. Snow pillows have measurement errors, logistical and transport problems for their installation, and can only measure an area of about 10 m 2 [15]. There are also other methods for measuring snow density, including proximal remote sensing (such as the GMON (GammaMONitor) snow water equivalent probe (Campbell Scientific Canada, Edmonton, AB, Canada)) [17], spatial remote sensing (microwave remote sensing) [18]. However, these methods have some drawbacks; for example, they do not measure the density of each snow layer that makes up the vertical stratigraphy of the snowpack, but only the average density of the snowpack. Recently, optical sensor data have been used as an alternative to monitor snow cover over large areas and have led to improved monitoring and management of this water resource [19]. To optimize the modeling process of the physical properties of snow and to develop new efficient models, a high spectral resolution is essential [20].
Hyperspectral imaging technology is an innovative approach based on spectroscopic analyses. It is fast, non-invasive, and facilitates real-time measurements [21,22], which can be used in conjunction with traditional measurement methods. This technology has proven to be effective for field, laboratory, and industrial applications [23,24]. It provides detailed information about the physical and chemical components of a scanned sample due to its high spectral and temporal resolution [23,25,26]. It has already been demonstrated that the near-infrared (NIR) spectrum is sensitive to the physical parameters of snow [27][28][29][30]. In fact, snow granulometry is clearly visible in the NIR and the short waves of infrared regions (SWIR) [31,32]. Eppanapelli et al. [33] found that the spectral reflectance of snow in the NIR is inversely proportional to the liquid water content in the snow. In addition, the absorption of ice in the NIR is very high [34], so the effect of impurities such as mineral dust and soot is negligible beyond 1000 nm wavelength. The above findings highlight the potential of hyperspectral NIR data to gather information on the physical properties of snow for modeling purposes [35].
Indeed, several models and approaches designed to model snow density based on spectral data are now available, but none has yet achieved a satisfactory performance [36,37]. This is probably due to the fact that most models are based on the assumption that density measurements can be modeled using the same function. Even though the spectral reflectance of snow in the NIR depends on density, it is expressed by the size and shape of the grains (granulometry) and the liquid water content in the snow [38]. Consequently, the optical properties of the physical parameters of snow influence one another and create a non-bijective (surjective (In mathematics, a surjective function (also known as surjection, or onto function) is a function f that maps an element x to every element y; that is, for every y, there is an x such that f(x) = y.)) relationship between snow density and reflectance in the NIR, resulting in poor correlations.
Recently, it has been demonstrated that three optical classes of snow with different degrees of metamorphosis (weakly to moderately metamorphosed (WMM), moderately to highly metamorphosed (MHM), and highly to very highly metamorphosed (HVM)) can be identified and discriminated against without prior recognition, based only on NIR hyperspectral data [39]. This study showed that the spectra of snow density are similar within the same optical class and significantly different from one optical class to another. In other words, densities expressed in terms of grain size, shape, and spectral response were discriminated and grouped into three different homogeneous subclasses [39]. With this finding, it is possible to train estimators specific to the identified homogeneous classes, which are governed by a bijective (In mathematics, a bijection, bijective function, one-to-one correspondence, or invertible function, is a function between the elements of two sets, where each element of one set is paired with exactly one element of the other set, and each element of the other set is paired with exactly one element of the first set. There are no unpaired elements.) relation between the density and the hyperspectral NIR reflectance.
The objective of this study is to develop a hybrid model (HM) to estimate the snow density using proximal NIR hyperspectral data. The HM is a combination of a classifier and specific estimators associated with three density classes (WMM, MHM, and HVM). The HM was calibrated and validated using a data set collected from a sampling site located in Quebec City, Canada. The performance of the HM was assessed using the leave-oneout cross-validation technique and independent validation data using a systematic data splitting technique. Four statistical evaluation indices (the coefficient of determination (R 2 ), root mean square error (RMSE), the bias (BIAS), and Nash-criterion (NASH)) were used to assess the model's performance.

Study Area
An extensive field campaign was conducted during the winter of (19 January-27 March) 2018, (10 January-3 April) 2019, and (29 January-10 March) 2020. The sampling zone is located in Quebec City (Canada) on the premises of INRS's (National Institute of Scientific Research) technology park (46 • 47 43.22" north latitude and −71 • 18 10" west longitude; Figure 1). An open area of approximately 20 m × 5 m located in a hardwood forest was selected as a sampling site. Measurements were collected in the morning between 8:00 a.m. and 12:00 p.m. on sunny and windless days. The climate in Quebec City is characterized by winter temperatures ranging between −10 • C and −25 • C, with significant snow accumulation. The snowpack in this region is dry from January to mid-February due to low temperatures and becomes wet in March as air temperature, day length, and radiation increase. to-one correspondence, or invertible function, is a function between the elements of two sets, where each element of one set is paired with exactly one element of the other set, and each element of the other set is paired with exactly one element of the first set. There are no unpaired elements.) relation between the density and the hyperspectral NIR reflectance. The objective of this study is to develop a hybrid model (HM) to estimate the snow density using proximal NIR hyperspectral data. The HM is a combination of a classifier and specific estimators associated with three density classes (WMM, MHM, and HVM). The HM was calibrated and validated using a data set collected from a sampling site located in Quebec City, Canada. The performance of the HM was assessed using the leaveone-out cross-validation technique and independent validation data using a systematic data splitting technique. Four statistical evaluation indices (the coefficient of determination (R 2 ), root mean square error (RMSE), the bias (BIAS), and Nash-criterion (NASH)) were used to assess the model's performance.

Study Area
An extensive field campaign was conducted during the winter of 2018 (19 January-27 March), 2019 (10 January-3 April), and 2020 (29 January-10 March). The sampling zone is located in Quebec City (Canada) on the premises of INRS's (National Institute of Scientific Research) technology park (46°47′43.22″ north latitude and −71°18′10″ west longitude; Figure 1). An open area of approximately 20 m × 5 m located in a hardwood forest was selected as a sampling site. Measurements were collected in the morning between 8:00 a.m. and 12:00 p.m. on sunny and windless days. The climate in Quebec City is characterized by winter temperatures ranging between −10 °C and −25 °C, with significant snow accumulation. The snowpack in this region is dry from January to mid-February due to low temperatures and becomes wet in March as air temperature, day length, and radiation increase.

In-Situ Data Collection
Two types of instruments were used to generate the calibration and validation database. Snow samples were collected using a rectangular core sampler designed and built by INRS's remote sensing team (Figure 2), and optical properties were measured using a proximal acquisition station ( Figure 3). The data acquired with these two instruments were collected at the same place and time.

In-Situ Data Collection
Two types of instruments were used to generate the calibration and validation database. Snow samples were collected using a rectangular core sampler designed and built by INRS's remote sensing team (Figure 2), and optical properties were measured using a proximal acquisition station ( Figure 3). The data acquired with these two instruments were collected at the same place and time.

In-Situ Data Collection
Two types of instruments were used to generate the calibration and validation database. Snow samples were collected using a rectangular core sampler designed and built by INRS's remote sensing team (Figure 2), and optical properties were measured using a proximal acquisition station ( Figure 3). The data acquired with these two instruments were collected at the same place and time.   The rectangular core sampler (or corer) was used to extract the vertical stratigraphy of the snowpack (Figure 4a). The corer is composed of a metallic inner component and a plastic cover with a triangular sawtooth cutting part at the end. This design allows the extraction of the vertical profile of the snowpack with all its metamorphosed variants and no loss of snow ( Figure 4b). The corer is graduated to measure the height (cm) of the profile and to differentiate the homogeneous layers composing the vertical profile. The rectangular core sampler (or corer) was used to extract the vertical stratigraphy of the snowpack (Figure 4a). The corer is composed of a metallic inner component and a plastic cover with a triangular sawtooth cutting part at the end. This design allows the extraction of the vertical profile of the snowpack with all its metamorphosed variants and no loss of snow ( Figure 4b). The corer is graduated to measure the height (cm) of the profile and to differentiate the homogeneous layers composing the vertical profile.  The proximal acquisition station used in this work is equipped with a hyperspectral camera (PIKA NIR from RESONON Company) boarded on a linear translation plate which allows for fast acquisition of spectral data ( Figure 3). The camera measures the NIR spectral reflectance for wavelengths between 900 nm and 1700 nm with a spectral resolution of 5.5 nm for 148 spectral bands. The station also contains a halogen lighting system, a mounting tower, a mobile platform, a data acquisition software (Spectronon Pro by Resonon Inc., Bozeman, MT, USA), and proximal acquisition lenses ( Figure 3).
The snow sample is collected by pushing the corer through the surface of the snowpack until the serrated cutting end reaches the ground. This ensures the recovery of both surface snow and snow from deeper layers (Figure 4a). Once the snow is extracted ( Figure  4b), it is scanned with the NIR hyperspectral camera (Figure 5a), and the generated image is analyzed by the Spectronon Pro software (Figure 5b) [39]. The latter identifies the homogeneous layers previously measured in the field and analyzes the spectral responses of each one. After this, each homogeneous layer is carefully removed from the corer, and its physical properties are characterized. The layers are then classified according to the International classification of seasonal ground snow [40]. To do so, the size and type of snow grains are identified using a millimeter grid map and observed with a magnifying glass (10×). Finally, the isolated layers are weighed, and the density of each one is calculated [39]. All observations and measurements were performed by the same person to ensure consistency. The proximal acquisition station used in this work is equipped with a hyperspectral camera (PIKA NIR from RESONON Company) boarded on a linear translation plate which allows for fast acquisition of spectral data ( Figure 3). The camera measures the NIR spectral reflectance for wavelengths between 900 nm and 1700 nm with a spectral resolution of 5.5 nm for 148 spectral bands. The station also contains a halogen lighting system, a mounting tower, a mobile platform, a data acquisition software (Spectronon Pro by Resonon Inc., Bozeman, MT, USA), and proximal acquisition lenses ( Figure 3).
The snow sample is collected by pushing the corer through the surface of the snowpack until the serrated cutting end reaches the ground. This ensures the recovery of both surface snow and snow from deeper layers (Figure 4a). Once the snow is extracted (Figure 4b), it is scanned with the NIR hyperspectral camera (Figure 5a), and the generated image is analyzed by the Spectronon Pro software (Figure 5b) [39]. The latter identifies the homogeneous layers previously measured in the field and analyzes the spectral responses of each one. After this, each homogeneous layer is carefully removed from the corer, and its physical properties are characterized. The layers are then classified according to the International classification of seasonal ground snow [40]. To do so, the size and type of snow grains are identified using a millimeter grid map and observed with a magnifying glass (10×). Finally, the isolated layers are weighed, and the density of each one is calculated [39]. All observations and measurements were performed by the same person to ensure consistency.

Methodological Approach
The methodological approach adopted is based on three main phases that are summarized in the flow chart illustrated in Figure 6:  The calibration of the classifier and the specific estimation of the HM The calibration was carried out in two steps: (1) Development of a machine learningbased algorithms using CART (classification and regression tree) method to discriminate between three different snow classes (WMM, MHM, and HVM) based on the NIR spectral data; The CART algorithm is a tree-based decision process. It uses a training sample set composed of a historical data set with pre-assigned classes for all observations (density in our case) and splitting variables (spectral bands in our case). These decision trees are spectrally based thresholds that are then used to classify new data [41]. And (2) Subdivision of the calibration database, based on the previously developed classifier, into three subdata sets. These sub-data sets were then used to train three specific estimators. For each specific estimator, all possible band ratios were first calculated and correlated with the snow density values. The correlation coefficients were then stored in a coefficient correlation matrix. The most highly correlated band ratios (R 2 > 0.5) were then integrated into a stepwise algorithm to train the calibration function. This exercise was also applied to the band differences and the normalized differences independently. Whether based on ratios, differences, or normalized differences, the spectral index most highly correlated with each specific estimator was retained at the end of the calibration step.  Evaluation of the specific estimators using the leave-one-out cross-validation (LOOCV) algorithm [42]; The specific estimators were assessed using the LOOCV method. This method consists of temporarily removing an entry from the database and using the rest of the database as calibration data, and then estimating the density of the removed measurement. This operation is repeated for the entire database. The main objective of this step was to calculate the relative bias of each specific estimator.

Methodological Approach
The methodological approach adopted is based on three main phases that are summarized in the flow chart illustrated in Figure 6:

Methodological Approach
The methodological approach adopted is based on three main phases that are summarized in the flow chart illustrated in Figure 6:  The calibration of the classifier and the specific estimation of the HM The calibration was carried out in two steps: (1) Development of a machine learningbased algorithms using CART (classification and regression tree) method to discriminate between three different snow classes (WMM, MHM, and HVM) based on the NIR spectral data; The CART algorithm is a tree-based decision process. It uses a training sample set composed of a historical data set with pre-assigned classes for all observations (density in our case) and splitting variables (spectral bands in our case). These decision trees are spectrally based thresholds that are then used to classify new data [41]. And (2) Subdivision of the calibration database, based on the previously developed classifier, into three subdata sets. These sub-data sets were then used to train three specific estimators. For each specific estimator, all possible band ratios were first calculated and correlated with the snow density values. The correlation coefficients were then stored in a coefficient correlation matrix. The most highly correlated band ratios (R 2 > 0.5) were then integrated into a stepwise algorithm to train the calibration function. This exercise was also applied to the band differences and the normalized differences independently. Whether based on ratios, differences, or normalized differences, the spectral index most highly correlated with each specific estimator was retained at the end of the calibration step.  Evaluation of the specific estimators using the leave-one-out cross-validation (LOOCV) algorithm [42]; The specific estimators were assessed using the LOOCV method. This method consists of temporarily removing an entry from the database and using the rest of the database as calibration data, and then estimating the density of the removed measurement. This operation is repeated for the entire database. The main objective of this step was to calculate the relative bias of each specific estimator.
The calibration of the classifier and the specific estimation of the HM The calibration was carried out in two steps: (1) Development of a machine learningbased algorithms using CART (classification and regression tree) method to discriminate between three different snow classes (WMM, MHM, and HVM) based on the NIR spectral data; The CART algorithm is a tree-based decision process. It uses a training sample set composed of a historical data set with pre-assigned classes for all observations (density in our case) and splitting variables (spectral bands in our case). These decision trees are spectrally based thresholds that are then used to classify new data [41]. And (2) Subdivision of the calibration database, based on the previously developed classifier, into three subdata sets. These sub-data sets were then used to train three specific estimators. For each specific estimator, all possible band ratios were first calculated and correlated with the snow density values. The correlation coefficients were then stored in a coefficient correlation matrix. The most highly correlated band ratios (R 2 > 0.5) were then integrated into a stepwise algorithm to train the calibration function. This exercise was also applied to the band differences and the normalized differences independently. Whether based on ratios, differences, or normalized differences, the spectral index most highly correlated with each specific estimator was retained at the end of the calibration step. Sens. 2021, 13, x FOR PEER REVIEW 8 of 20 Figure 6. Flowchart of the methodological approach. HM refers to hybrid model, S.E refers to specific estimator, and LOOCV refers to leave-one-out cross-validation.

Descriptive Analysis of Data
Domine, Taillandier, and Simpson [36] proposed snow grain size parameterizations based on snow density and grain shape, grouping snow samples into three main types (fresh, recent, and old). Based on this principle, we treated the snowpack as a succession of layers. Each layer is characterized by its own physical properties depending on its history and maturation stage of transformation [44]. The layers are reported in Table 1 and are grouped by grain size, grain type, and layer density, as defined by Pahaut [45]. Based on these physical properties, each layer of snow has been assigned to one of three general predefined snow classes (WMM, MHM, and HVM) following the classification of Pahaut [45] and the International classification of seasonal ground snow [40] protocols. For inventory purposes, a description and pictures of snow grains with their corresponding dates were taken for each layer. These were used later as real in-situ data when verification was needed.
A database of 114 layers was established (Table 1)

Methodological Approach
The methodological approach adopted is based on three main phases that are summarized in the flow chart illustrated in Figure 6:  The calibration of the classifier and the specific estimation of the HM The calibration was carried out in two steps: (1) Development of a machine learningbased algorithms using CART (classification and regression tree) method to discriminate between three different snow classes (WMM, MHM, and HVM) based on the NIR spectral data; The CART algorithm is a tree-based decision process. It uses a training sample set composed of a historical data set with pre-assigned classes for all observations (density in our case) and splitting variables (spectral bands in our case). These decision trees are spectrally based thresholds that are then used to classify new data [41]. And (2) Subdivision of the calibration database, based on the previously developed classifier, into three subdata sets. These sub-data sets were then used to train three specific estimators. For each specific estimator, all possible band ratios were first calculated and correlated with the snow density values. The correlation coefficients were then stored in a coefficient correlation matrix. The most highly correlated band ratios (R 2 > 0.5) were then integrated into a stepwise algorithm to train the calibration function. This exercise was also applied to the band differences and the normalized differences independently. Whether based on ratios, differences, or normalized differences, the spectral index most highly correlated with each specific estimator was retained at the end of the calibration step.  Evaluation of the specific estimators using the leave-one-out cross-validation (LOOCV) algorithm [42]; The specific estimators were assessed using the LOOCV method. This method consists of temporarily removing an entry from the database and using the rest of the database as calibration data, and then estimating the density of the removed measurement. This operation is repeated for the entire database. The main objective of this step was to Evaluation of the specific estimators using the leave-one-out cross-validation (LOOCV) algorithm [42]; The specific estimators were assessed using the LOOCV method. This method consists of temporarily removing an entry from the database and using the rest of the database as calibration data, and then estimating the density of the removed measurement. This operation is repeated for the entire database. The main objective of this step was to calculate the relative bias of each specific estimator.

Methodological Approach
The methodological approach adopted is based on three main phases that are summarized in the flow chart illustrated in Figure 6:  The calibration of the classifier and the specific estimation of the HM The calibration was carried out in two steps: (1) Development of a machine learningbased algorithms using CART (classification and regression tree) method to discriminate between three different snow classes (WMM, MHM, and HVM) based on the NIR spectral data; The CART algorithm is a tree-based decision process. It uses a training sample set composed of a historical data set with pre-assigned classes for all observations (density in our case) and splitting variables (spectral bands in our case). These decision trees are spectrally based thresholds that are then used to classify new data [41]. And (2) Subdivision of the calibration database, based on the previously developed classifier, into three subdata sets. These sub-data sets were then used to train three specific estimators. For each specific estimator, all possible band ratios were first calculated and correlated with the snow density values. The correlation coefficients were then stored in a coefficient correlation matrix. The most highly correlated band ratios (R 2 > 0.5) were then integrated into a stepwise algorithm to train the calibration function. This exercise was also applied to the band differences and the normalized differences independently. Whether based on ratios, differences, or normalized differences, the spectral index most highly correlated with each HM evaluation using independent data selected using the systematic split validation (SSV) technique.
Firstly, 25% of the data was removed from the initial database prior to the calibration of the MH for independent validation purposes. The validation data were systematically selected (the 4th data entries starting with low to high-density values) and set aside. This technique is known as systematic split validation (SSV). The remaining data (calibration database) was used to train the HM. Finally, to assess the overall accuracy of the HM, the independent dataset (the 25% removed by the SSV technique) was used. It is important to note, however, that it was necessary to go through both calibration steps to accurately assess the MH. Each sample was first assigned to one of three snow classes using the HM classifier. The snow density was then estimated using the specific estimator corresponding to the pre-assigned class.
Four statistical indices were used to assess both the specific estimators and the HM, which are presented in Equations (1)- (4). These are the coefficient of determination (R 2 ), defined as the squared value of the correlation coefficient, the root means square error (RMSE), the bias (BIAS), and the Nash-criterion (NASH). The NASH criterion [43] compares the values estimated by the model to the mean in-situ measurements, yielding a result between −∞ and 1.0 (inclusive). A negative NASH result means that it would be preferable to use the mean in-situ measurements rather than the model estimates, whereas values between 0.0 and 0.6 are generally considered acceptable levels of model performance. Model performance is satisfactory for values above 0.8, and the model is perfect when NASH = 1.0. The mathematical equations for these indices are as follows: where: n is the size of the dataset, M and Es are the measured and estimated density values, M and Es are the means of the measured and estimated values.

Descriptive Analysis of Data
Domine, Taillandier, and Simpson [36] proposed snow grain size parameterizations based on snow density and grain shape, grouping snow samples into three main types (fresh, recent, and old). Based on this principle, we treated the snowpack as a succession of layers. Each layer is characterized by its own physical properties depending on its history and maturation stage of transformation [44]. The layers are reported in Table 1 and are grouped by grain size, grain type, and layer density, as defined by Pahaut [45]. Based on these physical properties, each layer of snow has been assigned to one of three general predefined snow classes (WMM, MHM, and HVM) following the classification of Pahaut [45] and the International classification of seasonal ground snow [40] protocols. For inventory purposes, a description and pictures of snow grains with their corresponding dates were taken for each layer. These were used later as real in-situ data when verification was needed.
A database of 114 layers was established (Table 1) Table 1, we consider newly deposited snow on the snowpack to be characterized by small grain size and low density, while moderately and heavily metamorphosed snow was characterized by moderate and large grain size and density, respectively. The smaller number of layers in the WMM class is due to the difficulty of isolating this type of layer because weakly metamorphosed snow is very light and often thin.  Figure 7 shows the measured NIR reflectance spectrum of the 114 snow samples (layers) extracted from twenty-four snow cores (all field campaigns included). The shape of the spectral reflectance changes as a function of the physical properties of the snow (density, grain size, and grain shape), which depends on the state of transformation by the metamorphosis of each layer. In fact, each layer is characterized by its own spectral reflectance, which is characterized by the degree of reflection, absorption, and transmission in each wavelength of the NIR spectrum. This results in similarities between the reflectance spectra of layers with similar chemical composition and physical characteristics [46]. This is translated by a decrease of reflectance as the snow layer ages; WMM, MHM, and HVM are colored dark blue, sky blue, and pink, respectively. These findings show that it is possible to discriminate between the three predefined classes based only on the reflectance spectra of the layers. grain size and density, respectively. The smaller number of layers in the WMM class is due to the difficulty of isolating this type of layer because weakly metamorphosed snow is very light and often thin. Table 1. Distribution of snow density as a function of type and size of snow grains, for WMM, MHM, and HVM, respectively, (field data: 2018, 2019 and 2020). Precipitation particles (+), decomposing and fragmented precipitation particles (λ), rounded grains (•), faceted crystals (□), depth hoar (˄), melt forms (ᴼ). The different shades of gray show the levels of metamorphosis (low, moderate, high) of each recuperated snow layer.

Grain Size (mm)
Grain Type  Figure 7 shows the measured NIR reflectance spectrum of the 114 snow samples (layers) extracted from twenty-four snow cores (all field campaigns included). The shape of the spectral reflectance changes as a function of the physical properties of the snow (density, grain size, and grain shape), which depends on the state of transformation by the metamorphosis of each layer. In fact, each layer is characterized by its own spectral reflectance, which is characterized by the degree of reflection, absorption, and transmission in each wavelength of the NIR spectrum. This results in similarities between the reflectance spectra of layers with similar chemical composition and physical characteristics [46]. This is translated by a decrease of reflectance as the snow layer ages; WMM, MHM, and HVM are colored dark blue, sky blue, and pink, respectively. These findings show that it is possible to discriminate between the three predefined classes based only on the reflectance spectra of the layers.

Statistical Dependence between Density and Spectral Reflectance
Previous works have failed to establish a strong correlation between snow density and reflectance, neither in the visible [37,47] nor in the NIR [36,37,47] parts of the spectrum. The same poor correlations (maximum correlation of R 2 = 0.53) were also found in this study when using single bands, as shown by the graphical correlogram shown in Figure 8a. The R 2 is also presented as a correlation matrix using three spectral indices: normalized band differences (Figure 8b), band ratios (Figure 8c), and band differences (Figure 8d). The accuracy was relatively increased, and the results are consistent with other published works [33,34]. The best correlation for the normalized band differences, band ratios, and band differences was found with the spectral ranges of 1034-1051 nm (R 2 = 0.64), 1205-1309 nm (R 2 = 0.63), and 1644-1655 nm (R 2 = 0.56), respectively. The above results affirmed that the correlation between the snow density and the spectral indices is relatively moderate (R 2 =< 0.64).
trum. The same poor correlations (maximum correlation of R 2 = 0.53) were also found in this study when using single bands, as shown by the graphical correlogram shown in Figure 8a. The R 2 is also presented as a correlation matrix using three spectral indices: normalized band differences (Figure 8b), band ratios (Figure 8c), and band differences ( Figure  8d). The accuracy was relatively increased, and the results are consistent with other published works [33,34]. The best correlation for the normalized band differences, band ratios, and band differences was found with the spectral ranges of 1034-1051 nm (R 2 = 0.64), 1205-1309 nm (R 2 = 0.63), and 1644-1655 nm (R 2 = 0.56), respectively. The above results affirmed that the correlation between the snow density and the spectral indices is relatively moderate (R 2 =< 0.64). Based on the above results, it was decided to go through the calibration and validation steps with the highest correlation, which was the normalized band differences. It is Based on the above results, it was decided to go through the calibration and validation steps with the highest correlation, which was the normalized band differences. It is important to note, however, that these results are based on 75% of the initial database, as the remaining 25% was utilized for validation purposes. In order to avoid autocorrelation between independent variables, it was decided to correlate only the bands with the best associations (here called neo-variables for simplicity) with snow density. The result is shown in Figure 9a. The density is sensitive to wavelengths 974 nm and 990 nm, with an R 2 of 0.64. When the model was validated using the independent data set (25% removed), all statistical indices showed moderate accuracy in estimating snow density, particularly R 2 and Nash (R 2 = Nash = 0.66; Figure 9b). the remaining 25% was utilized for validation purposes. In order to avoid autocorrelation between independent variables, it was decided to correlate only the bands with the best associations (here called neo-variables for simplicity) with snow density. The result is shown in Figure 9a. The density is sensitive to wavelengths 974 nm and 990 nm, with an R 2 of 0.64. When the model was validated using the independent data set (25% removed), all statistical indices showed moderate accuracy in estimating snow density, particularly R 2 and Nash (R 2 = Nash = 0.66; Figure 9b). Figure 9. Results of (a) calibration using normalized band differences (75% of the database), and (b) validation (25% of the database).
The aforementioned analysis confirms that the reciprocal correlation between the spectral reflectance in the NIR and snow density is probably the result of other changes in the aging process, such as temperature, day length, radiation, etc. [35]. Indeed, despite the strong correlation between the decrease in reflectance and the increase in density ( Figure  7), it is difficult to evaluate the extent to which the decrease should be attributed to the underlying physical cause of densification, and this is where the difficulty of experimentally relating spectral reflectance (using single bands or based on neo-variables) to snow density lies.

HM Classifier Calibration
Although several density models exist, most assume that snow density can be modeled using a single function (such as the one developed above). The use of a single function, even if it is multi-varied, is still confronted with the inherent surjective relationship between snow density and spectral reflectance. Therefore, replacing the surjective relation with several bijective relations could be the answer to increase the accuracy of snow density modeling. This could be achieved by establishing a strong correlation between the decrease in reflectance and the increase in density. To validate this hypothesis, a classification algorithm was applied to the calibration database to discriminate between the three predefined snow classes (WMM, MHM, and HVM) based on spectral reflectance and, hopefully, to increase the relationship between snow density and spectral reflectance in each class; converting this surjective relationship into three bijective relationships.
The results of in-situ measurement classification, presented in Table 1, were used to develop the HM's classifier using the CART algorithm. CART results are represented in the form of an intuitive illustration ( Figure 10). Two wavelengths (1024 nm and 1161 nm) Figure 9. Results of (a) calibration using normalized band differences (75% of the database), and (b) validation (25% of the database).
The aforementioned analysis confirms that the reciprocal correlation between the spectral reflectance in the NIR and snow density is probably the result of other changes in the aging process, such as temperature, day length, radiation, etc. [35]. Indeed, despite the strong correlation between the decrease in reflectance and the increase in density (Figure 7), it is difficult to evaluate the extent to which the decrease should be attributed to the underlying physical cause of densification, and this is where the difficulty of experimentally relating spectral reflectance (using single bands or based on neo-variables) to snow density lies.

HM Classifier Calibration
Although several density models exist, most assume that snow density can be modeled using a single function (such as the one developed above). The use of a single function, even if it is multi-varied, is still confronted with the inherent surjective relationship between snow density and spectral reflectance. Therefore, replacing the surjective relation with several bijective relations could be the answer to increase the accuracy of snow density modeling. This could be achieved by establishing a strong correlation between the decrease in reflectance and the increase in density. To validate this hypothesis, a classification algorithm was applied to the calibration database to discriminate between the three predefined snow classes (WMM, MHM, and HVM) based on spectral reflectance and, hopefully, to increase the relationship between snow density and spectral reflectance in each class; converting this surjective relationship into three bijective relationships.
The results of in-situ measurement classification, presented in Table 1, were used to develop the HM's classifier using the CART algorithm. CART results are represented in the form of an intuitive illustration ( Figure 10). Two wavelengths (1024 nm and 1161 nm) were selected as the best classifier thresholds, which allowed us to calibrate the classifier of the HM. Dozier and Painter [35] demonstrated that the spectral reflectance of snow displays an ice absorption feature centered at 1030 nm that can be used to estimate snow grain size from hyperspectral data. Roy et al. [48] used this ice absorption feature to distinguish natural snow from compacted or metamorphosed snow. The wavelengths identified by the CART are coherent with the works cited above. According to Figure 10, the band centered at 1024 nm allows to differentiate between the highly metamorphosed snow class (composed of snow grains of size > 2 mm) and the other two classes of weakly and moderately metamorphosed snow and the band centered at 1161 nm allows to discriminate between the moderately and weakly metamorphosed snow (composed of snow grains of size < 1 mm). tified by the CART are coherent with the works cited above. According to Figure 10, the band centered at 1024 nm allows to differentiate between the highly metamorphosed snow class (composed of snow grains of size > 2 mm) and the other two classes of weakly and moderately metamorphosed snow and the band centered at 1161 nm allows to discriminate between the moderately and weakly metamorphosed snow (composed of snow grains of size < 1 mm).

Calibration of the Specific Models
The calibrated classifier was then used to subdivide the calibration database into three subclasses based on spectra. Fifteen samples were used as a training dataset for the WMM class, 44 for the MHM class, and 27 for the HVM class. Each subclass dataset was thereafter used to calibrate its own specific estimator. The calibration step was performed using multivariate stepwise regression in which the independent variables were selected using an automatic procedure that typically followed a sequence of F-tests. Stepwise regression was used in a direct selection mode, starting with none of the variables in the model, testing them one by one, and including only statistically significant variables [49]. For comparison purposes, the calibration and validation steps were performed exactly as indicated in the section above. However, only significant results are presented to alleviate the text. As expected, the correlation between snow density and spectral reflectance was higher when each class was considered separately. Snow density showed a wider range of sensitivity to neo-variables, which was different for band differences and normalized band differences, depending on the class ( Figure 11). These findings support the assumption that destroying the surjective relationship could lead to the creation of a bijective relationship in each class and consequently enhancing the model's accuracy. Also, the highest correlations are not necessarily centered at the same wavelength for each class. This is a demonstration that each class has its own response to the spectral reflectance and the optimal wavelengths to detect snow density shift depending on the class. This is to be expected since the density depends on the grain size, the grain type, and the liquid water

Calibration of the Specific Models
The calibrated classifier was then used to subdivide the calibration database into three subclasses based on spectra. Fifteen samples were used as a training dataset for the WMM class, 44 for the MHM class, and 27 for the HVM class. Each subclass dataset was thereafter used to calibrate its own specific estimator. The calibration step was performed using multivariate stepwise regression in which the independent variables were selected using an automatic procedure that typically followed a sequence of F-tests. Stepwise regression was used in a direct selection mode, starting with none of the variables in the model, testing them one by one, and including only statistically significant variables [49]. For comparison purposes, the calibration and validation steps were performed exactly as indicated in the section above. However, only significant results are presented to alleviate the text. As expected, the correlation between snow density and spectral reflectance was higher when each class was considered separately. Snow density showed a wider range of sensitivity to neo-variables, which was different for band differences and normalized band differences, depending on the class ( Figure 11). These findings support the assumption that destroying the surjective relationship could lead to the creation of a bijective relationship in each class and consequently enhancing the model's accuracy. Also, the highest correlations are not necessarily centered at the same wavelength for each class. This is a demonstration that each class has its own response to the spectral reflectance and the optimal wavelengths to detect snow density shift depending on the class. This is to be expected since the density depends on the grain size, the grain type, and the liquid water content. In addition, it is known that water strongly absorbs in the NIR [34]. Hence, the higher the water content of snow, the smaller the returned signal (and vice versa). content. In addition, it is known that water strongly absorbs in the NIR [34]. Hence, the higher the water content of snow, the smaller the returned signal (and vice versa). Results of the stepwise algorithm are presented in Figure 12. The three snow classes were linearly correlated, and the calibration functions were univariate ( Table 2). The neovariable retained for the WMM class was the difference between 1265 nm and 941 nm wavelengths (referred to as the spectral index of subtraction; SISUB) with an R 2 of 0.91. The neo-variable retained for the MHM class was the normalized difference between 1617 nm and 941 nm wavelengths (referred here as the spectral index of normalized difference; SINOR) with an R 2 of 0.80. The neo-variable retained for the HVM class was also the difference between 1424 nm and 1188 nm wavelength (SISUB) with an R 2 of 0.84. Results of the stepwise algorithm are presented in Figure 12. The three snow classes were linearly correlated, and the calibration functions were univariate ( Table 2). The neovariable retained for the WMM class was the difference between 1265 nm and 941 nm wavelengths (referred to as the spectral index of subtraction; SI SUB ) with an R 2 of 0.91. The neo-variable retained for the MHM class was the normalized difference between 1617 nm and 941 nm wavelengths (referred here as the spectral index of normalized difference; SI NOR ) with an R 2 of 0.80. The neo-variable retained for the HVM class was also the difference between 1424 nm and 1188 nm wavelength (SI SUB ) with an R 2 of 0.84. Remote Sens. 2021, 13, x FOR PEER REVIEW 14 of 20 These results are coherent with other works. For example, Negi et al. [46] found that for changes in liquid water content, grain size, and metamorphosis level (aging) of snow, the maximum spectral variations are observed around the shorter wavelengths of the NIR spectrum. Gallet et al. [50] used the NIR spectrum to determine the size of snow grains. They found that the band centered at 1310 nm, which corresponds to the central part of the NIR spectrum, is sensitive to small snow grains, which are usually of low to medium density, and that the band centered at 1550 nm, which corresponds to the higher part of the NIR spectrum, is sensitive to large snow grains, where the snow is generally denser. The strong correlation between in situ density measurements and spectral indices (R 2 > 0.80) for all three specific estimators illustrates the importance of dividing the input data into several classes of snow (three in this study), based on a certain degree of aging. For all specific estimators, the estimation of snow density was limited to short (941 nm, 1265 nm, and 1188 nm) and long (1424 nm and 1617 nm) NIR wavelengths, which are the spectral regions most sensitive to the physical properties of snow. The mathematical expressions of the three specific estimators are summarized in Table 2. These results are coherent with other works. For example, Negi et al. [46] found that for changes in liquid water content, grain size, and metamorphosis level (aging) of snow, the maximum spectral variations are observed around the shorter wavelengths of the NIR spectrum. Gallet et al. [50] used the NIR spectrum to determine the size of snow grains. They found that the band centered at 1310 nm, which corresponds to the central part of the NIR spectrum, is sensitive to small snow grains, which are usually of low to medium density, and that the band centered at 1550 nm, which corresponds to the higher part of the NIR spectrum, is sensitive to large snow grains, where the snow is generally denser. The strong correlation between in situ density measurements and spectral indices (R 2 > 0.80) for all three specific estimators illustrates the importance of dividing the input data into several classes of snow (three in this study), based on a certain degree of aging. For all specific estimators, the estimation of snow density was limited to short (941 nm, 1265 nm, and 1188 nm) and long (1424 nm and 1617 nm) NIR wavelengths, which are the spectral regions most sensitive to the physical properties of snow. The mathematical expressions of the three specific estimators are summarized in Table 2. The LOOCV is a technique that consists of temporarily removing certain entries from a database and using the rest as calibration data, and then estimating the value of the removed entry. This operation is repeated for the entire database, resulting in an estimate for the entire database, hence allowing the comparison between estimated and measured values (Equations (1)-(4)). Figure 13 shows  The LOOCV is a technique that consists of temporarily removing certain entries from a database and using the rest as calibration data, and then estimating the value of the removed entry. This operation is repeated for the entire database, resulting in an estimate for the entire database, hence allowing the comparison between estimated and measured values (Equations (1)-(4)). Figure 13 shows  The LOOCV was not intended to evaluate the HM's snow density estimates, and its main objective was to quantify the BIAS of each specific estimator. In fact, each modeling process is tainted by two sources of error: (1) a random error quantified by the RMSE, and (2) a systematic error quantified by the BIAS. The first one is not correctable, but the second one is adjustable and can be removed/added to the estimates during the modeling process. This part of the study was aimed at quantifying the systematic error and using it to adjust the final model.

Evaluation of the Hybrid Model Using SSV
The HM was evaluated using independent data. However, it is important to note that in order to use SSV data for validation, it was necessary to go through all the calibration steps:

Methodological Approach
The methodological approach adopted is based on three main phases that are summarized in the flow chart illustrated in Figure 6:  The calibration of the classifier and the specific estimation of the HM The calibration was carried out in two steps: (1) Development of a machine learningbased algorithms using CART (classification and regression tree) method to discriminate between three different snow classes (WMM, MHM, and HVM) based on the NIR spectral data; The CART algorithm is a tree-based decision process. It uses a training sample set composed of a historical data set with pre-assigned classes for all observations (density in our case) and splitting variables (spectral bands in our case). These decision trees are spectrally based thresholds that are then used to classify new data [41]. And (2) Subdivision of the calibration database, based on the previously developed classifier, into three subdata sets. These sub-data sets were then used to train three specific estimators. For each specific estimator, all possible band ratios were first calculated and correlated with the snow density values. The correlation coefficients were then stored in a coefficient correlation matrix. The most highly correlated band ratios (R 2 > 0.5) were then integrated into a stepwise algorithm to train the calibration function. This exercise was also applied to the band differences and the normalized differences independently. Whether based on ratios, differences, or normalized differences, the spectral index most highly correlated with each specific estimator was retained at the end of the calibration step.  Evaluation of the specific estimators using the leave-one-out cross-validation (LOOCV) algorithm [42]; The specific estimators were assessed using the LOOCV method. This method consists of temporarily removing an entry from the database and using the rest of the database as calibration data, and then estimating the density of the removed measurement. This operation is repeated for the entire database. The main objective of this step was to calculate the relative bias of each specific estimator.

Methodological Approach
The methodological approach adopted is based on three main phases that are summarized in the flow chart illustrated in Figure 6:  The calibration of the classifier and the specific estimation of the HM The calibration was carried out in two steps: (1) Development of a machine learningbased algorithms using CART (classification and regression tree) method to discriminate between three different snow classes (WMM, MHM, and HVM) based on the NIR spectral data; The CART algorithm is a tree-based decision process. It uses a training sample set composed of a historical data set with pre-assigned classes for all observations (density in our case) and splitting variables (spectral bands in our case). These decision trees are spectrally based thresholds that are then used to classify new data [41]. And (2) Subdivision of the calibration database, based on the previously developed classifier, into three subdata sets. These sub-data sets were then used to train three specific estimators. For each specific estimator, all possible band ratios were first calculated and correlated with the snow density values. The correlation coefficients were then stored in a coefficient correlation matrix. The most highly correlated band ratios (R 2 > 0.5) were then integrated into a stepwise algorithm to train the calibration function. This exercise was also applied to the band differences and the normalized differences independently. Whether based on ratios, differences, or normalized differences, the spectral index most highly correlated with each specific estimator was retained at the end of the calibration step.  Evaluation of the specific estimators using the leave-one-out cross-validation (LOOCV) algorithm [42]; The specific estimators were assessed using the LOOCV method. This method consists of temporarily removing an entry from the database and using the rest of the database as calibration data, and then estimating the density of the removed measurement. This operation is repeated for the entire database. The main objective of this step was to calculate the relative bias of each specific estimator.
Estimate the density using the specific estimators corresponding to the preassigned class;

Methodological Approach
The methodological approach adopted is based on three main phases that are summarized in the flow chart illustrated in Figure 6:  The calibration of the classifier and the specific estimation of the HM The calibration was carried out in two steps: (1) Development of a machine learningbased algorithms using CART (classification and regression tree) method to discriminate between three different snow classes (WMM, MHM, and HVM) based on the NIR spectral data; The CART algorithm is a tree-based decision process. It uses a training sample set composed of a historical data set with pre-assigned classes for all observations (density in our case) and splitting variables (spectral bands in our case). These decision trees are spectrally based thresholds that are then used to classify new data [41]. And (2) Subdivision of the calibration database, based on the previously developed classifier, into three subdata sets. These sub-data sets were then used to train three specific estimators. For each specific estimator, all possible band ratios were first calculated and correlated with the snow density values. The correlation coefficients were then stored in a coefficient correlation matrix. The most highly correlated band ratios (R 2 > 0.5) were then integrated into a stepwise algorithm to train the calibration function. This exercise was also applied to the band differences and the normalized differences independently. Whether based on ratios, differences, or normalized differences, the spectral index most highly correlated with each specific estimator was retained at the end of the calibration step.  Evaluation of the specific estimators using the leave-one-out cross-validation (LOOCV) algorithm [42]; The specific estimators were assessed using the LOOCV method. This method consists of temporarily removing an entry from the database and using the rest of the database as calibration data, and then estimating the density of the removed measurement. This operation is repeated for the entire database. The main objective of this step was to calculate the relative bias of each specific estimator.
Correct the estimates' BIAS; Figure 14 shows the results of the evaluation of the HM's ability to estimate the snow density. The points are well distributed along the 1:1 line, from very low to high-density values. The value of the NASH index (R 2 = 0.93) highlights the good performance of the model. A satisfactory value of RMSE (31.48 kg m −3 ) is also achieved. As these results show, the performance of the HM is satisfactory, especially with a dataset covering a three-winter period. This period may seem long, but with global warming and climate change affecting climatic systems around the world, the seasonal snowpack is expected to also be changing from year to year. This is why finding an independent dataset to assess the performance of the HM was very challenging. Despite this, the evaluation shows that the HM accurately modeled the snow density using optical data. Contrary to previous studies claiming that there is no strong correlation between spectral reflectance and snow density, this study has shown that it is possible to estimate snow density if the right modeling strategy is found.
Remote Sens. 2021, 13, x FOR PEER REVIEW 17 of 20 Figure 14. Evaluation of the regional hybrid model using the SSV data.

Conclusions
The objective of this study was to test the performance of a hybrid model (HM) designed to estimate the density of the seasonal snowpack using hyperspectral NIR imaging (900-1700 nm) at a spectral resolution of 5.5 nm. The hybrid model is a combination of a classifier and three specific estimators (weakly to moderately metamorphosed snow (WMM), moderately to highly metamorphosed snow (MHM), and highly to very highly metamorphosed snow (HVM)). The hybrid model was evaluated at two levels: using the leave-one-out cross-validation (LOOCV) algorithm and using the systematic division validation technique (SSV). The LOOCV technique was used to assess the three specific estimators, and the SSV data were used to assess the performance of the HM.
The calibration step, based on a stepwise multivariate regression, showed that the Figure 14. Evaluation of the regional hybrid model using the SSV data.
The implementation of the approach developed in this work can constitute a real advance in the determination and monitoring of the various processes that lead to the continuous evolution and regular monitoring of the density of the seasonal snow cover. The results also highlight many of the difficulties inherent in measuring the density of snow on the ground, which presents considerable lateral variability (accumulation of layers of different densities). This large lateral variability is common to all snowpacks, and measuring only the average density of the latter does not represent the overall vertical profile. The approach described here, when fully matured, provides a means of obtaining the necessary data for these purposes.

Conclusions
The objective of this study was to test the performance of a hybrid model (HM) designed to estimate the density of the seasonal snowpack using hyperspectral NIR imaging (900-1700 nm) at a spectral resolution of 5.5 nm. The hybrid model is a combination of a classifier and three specific estimators (weakly to moderately metamorphosed snow (WMM), moderately to highly metamorphosed snow (MHM), and highly to very highly metamorphosed snow (HVM)). The hybrid model was evaluated at two levels: using the leave-one-out cross-validation (LOOCV) algorithm and using the systematic division validation technique (SSV). The LOOCV technique was used to assess the three specific estimators, and the SSV data were used to assess the performance of the HM.
The calibration step, based on a stepwise multivariate regression, showed that the three classes of snow are sensitive to different regions of the NIR spectrum, limited to the short and long wavelengths. The WMM was sensitive to the wavelengths at 1265 nm and 941 nm, the MHM was sensitive to wavelengths at 1617 nm and 941 nm, and the HVM was sensitive to wavelengths at 1424 nm and 1188 nm. The LOOCV technique highlighted that the specific estimators of all classes tend to slightly overestimate the snow density (BIAS < 0.1 kg·m −3 ). When the HM was challenged with SSV data, the modeling results were satisfactory with an R 2 = Nash = 0.93, and the snow density was slightly underestimated (BIAS = 1.03 kg·m −3 ).
The objective of this study was to develop a method based on the optical properties of snow to be used conjointly with conventional density measurement methods with the aim of alleviating field operations. The critical step in estimating snow density using the HM is the selection of the final specific estimator. Indeed, classification algorithms (such as CART) are known to be local and unstable. This instability can significantly affect the accuracy of the density using the specific estimators of the HM. In other words, for an ideal modeling process using the HM, the sample to be modeled must be well classified so that the specific estimator corresponding to that class is used for optimal density estimation. Otherwise, a wrong specific estimator will be selected, and consequently, the estimation will not be optimal, which could affect the accuracy. For example, for a measured density of 581 kg m −3 (classified as HVM), the relative errors vary by 5%, 39%, and 75% when estimated using specific estimators of HVM, MHM, and WMM, respectively. On the other hand, another obstacle associated with this method is the correct selection of homogeneous snow layers in the field and on the recovered hyperspectral image. For this reason, additional field campaigns need to be conducted to collect more data to overcome this weakness and allow proper field implementation. The HM provides an improved tool to monitor the evolution of seasonal snowpacks, with a satisfactory level of performance even for low to moderate snow densities. We conclude that our results are an important first step toward the development of an effective method for continuous monitoring of snow density profiles in the field.