What is the Trade-O ﬀ between Snowpack Stratiﬁcation and Simulated Snow Water Equivalent in a Physically-Based Snow Model?

: In Nordic watersheds, estimation of the dynamics of snow water equivalent (SWE) represents a major step toward a satisfactory modeling of the annual hydrograph. For a multilayer, physically-based snow model like MASiN (Mod è le Autonome de Simulation de la Neige), the number of modeled snow layers can a ﬀ ect the accuracy of the simulated SWE. The objective of this study was to identify the maximum number of snow layers (MNSL) that would deﬁne the trade-o ﬀ between snowpack stratiﬁcation and SWE modeling accuracy. Results indicated that decreasing the MNSL reduced the SWE modeling accuracy since the thermal energy balance and the mass balance were less accurately resolved by the model. Nevertheless, from a performance standpoint, SWE modeling can be accurate enough with a MNSL of two (2), with a substantial performance drop for a MNSL value of around nine (9). Additionally, the linear correlation between the values of the calibrated parameters and the MNSL indicated that reducing the latter in MASiN increased the fresh snow density and the settlement coe ﬃ cient, while the maximum radiation coe ﬃ cient decreased. In this case, MASiN favored the melting process, and thus the homogenization of snow layers occurred from the top layers of the snowpack in the modeling algorithm.


Introduction
In northern watersheds, snowfall constitutes a significant proportion of the total precipitation [1,2]. When rainfall happens at a low rate, water infiltrates until the soil becomes saturated, at which point surface runoff occurs. In contrast, when snowfall takes place, water is stored on the ground surface in solid form, modifying the soil water dynamic throughout the watershed [3]. The snowpack starts melting as it absorbs the amount of energy required for phase change to occur. Although this process can happen partially during winter, the atmospheric warming during spring will gradually melt the entire snowpack. Thus, accumulation and gradual snowmelt processes temporally alter a significant proportion of the total precipitation to flow toward the river network or recharge the groundwater [4][5][6].
Various snow modeling approaches exist. Some models such as CROCUS [7,8] emphasize the internal properties of the snow cover (e.g., snowpack stratigraphy), which can provide helpful information with respect to snowpack stability. This information becomes valuable, for instance, to predict snow avalanches in mountainous regions. Other types of snow models solely focus on snow water equivalent (SWE) in order to determine the timing and amount of melt. Snow models can be unique set of valid parameters at nearby snow sites. This is an advantage when it comes to pairing MASiN with a hydrological model for flow modeling, particularly when the hydrological model is semi-distributed or distributed because it is necessary to model the snow melt dynamics in different parts of a watershed.
Before pairing MASiN with any hydrological model such as HYDROTEL [25,26] for pragmatic reasons, it is noteworthy to account for the computational budget of model calibration and the operational application for inflow forecasting. Both procedures are indeed directly affected by the value of the maximum number of snow layers (MNSL). In MASiN, the MNSL value can be specified by the user. Mas et al. [24] considered a MNSL value of 70 to discretize the mass-energy transfer between snow layers. In HYDROTEL, the watershed is discretized into hillslopes referred to as relatively homogeneous hydrological units (RHHUs). Similarly, HYDROTEL discretizes the river network into different reaches that are fed by two or three RHHUs, the latter in the case of headwater reach (i.e., stream order 1 in a Strahler reference system). The number of RHHUs is not bounded, thus, from a computational standpoint, this would majorly impact the computational budget of MASiN for snowmelt simulation given the MNSL to simulate per RHHU.
Since MASiN is a physically-based snow model, the MNSL value can be analogous to the number of snow layers observed in nature. For example, in Antarctica, Arndt et al. [27] estimated the variability of snow cover properties at various sites. They observed that the plurennial snow cover had eight layers, while the seasonal snow cover had four layers. In Italy, Monti et al. [28] observed between seven and 20 layers at different times based on snow grain type, size, hardness, and density. Brun et al. [8] studied the operational forecasting of avalanches by simulating snow cover stratification and observed up to 13 layers with different types of snow grains in the French Alps. Finally, Armstrong [29] studied the compressive stress in the snow cover between neighboring layers and observed six different layers in Colorado, USA. Thus, in this study, in order to consider a coherent maximum number of layers to simulate, the MNSL in MASiN was set between one and 20 in comparison with the full-extent configuration of the original model, that is, a maximum of 70 snow layers.
Identifying the MNSL value is identical to modifying the total number of interactions in the model by reducing the number of computations of the energy and mass balances. In order to have an acceptable computational budget and to contribute to the development of this snow model, the flexibility, provided by this parameter, was investigated. This was done by evaluating the influence of the MNSL on model performance for snowmelt (expressed in snow water equivalent; SWE) modeling. When reducing the MNSL, bearing in mind a reasonable computational budget while maintaining a decent level of estimation accuracy, we assumed that this knowledge might prove to be helpful when pairing MASiN with any other distributed or semi-distributed hydrological model. To this end, parameter sets of MASiN were calibrated using the dynamically dimensioned search (DDS) optimization algorithm [30] at various snow stations by setting the MNSL value. It is also possible to compare the performance of MASiN, depending on the choice of the MNSL value, to simulate and compare its impact on the calibrated value of each parameter. Furthermore, modifying the MNSL value changes the number of interactions between snow layers, which in turn may impact on the value of another calibration parameter in the model. Hence, this study also investigated the extent to which other parameters were influenced when modifying the MNSL for SWE simulation. Based on this analysis, the user could thus adjust the MNSL value when it is required to amplify or reduce any modeled physical processes to improve SWE modeling performance. Moreover, while most modeling studies have generally focused on snowpack stratification from a snow properties point of view such as thermal conductivity or snow density (for instance, [17]) or on snow characteristics for preventing avalanche hazard (for instance, [31]), our study focused on snowpack stratification from a hydrological modeling point of view, and to our knowledge, this represents an original contribution.

Study Watersheds
Two watersheds were considered to test the model in different regions with diverse physical and regional characteristics. The locations of the watersheds within Canada are shown in Figure 1. Necopastic is a sub-watershed of the James Bay watershed with an area of 244 km 2 . The watershed is covered by coniferous forest (55%), bogs (35%), and outcrops and open water (10%). Data from the meteorological station located at the La Grande Rivière airport indicate that the average annual precipitation is 697 mm, 35% of which falls in the form of snow (true for 1981-2010 data from Environment and Climate Change Canada [32]). During this period, the average annual temperature over the watershed is −2.9 • C, and minimum and maximum average monthly temperatures are −28 • C in January and 20.4 • C in July, respectively. The Upper Yukon, on the other hand, has an area close to 20,000 km 2 , covered by coniferous forests (41.1%), bare ground/grass/shrub (42.4%), and ice and open water (16.4%). Data from the meteorological station located at the Whitehorse airport indicate that the average annual precipitation is 262 mm, 39% of which falls in the form of snowfall (true for 1981-2010 data from Environment and Climate Change Canada [33]. During this period, the average annual temperature over the watershed is −0.1 • C, and the minimum and maximum average monthly temperatures are −19.2 • C in January and 20.6 • C in July, respectively.

Meteorological Data
Input data including hourly air temperature, relative humidity, wind speed, and daily precipitation provided by the meteorological stations located near the snow stations are provided in Table 1 (see Figure 1 for the station locations). Ground precipitation measurements were intermittent during the observation period in the Upper Yukon, and, therefore, the precipitation time series had to be reconstructed for each snow station. For this, daily increases in observed SWE time series were considered as an input to the snowpack modeling. The source of these observed water inputs was assumed to be solely provided by the precipitation caused because of a lack of data about wind-induced snow drifting, and therefore, daily increase in the observed SWE time series were obtained from daily precipitation amounts of equal water depth.

Snow Data
Outputs from MASiN were compared to the observed SWE time series from snow stations, also called GMON (Gamma MONitor) stations. A GMON station is an automatic equipment developed by Choquette et al. [35] that continuously measures the snow cover (in mm of SWE) by estimating the natural ground gamma emission absorbed by the water content of the snow cover. As this study was based on the potential pairing of MASiN with HYDROTEL, the modeling performance was estimated using daily GMON time series for the different stations presented in Table 2 and shown in Figure 1. It is noteworthy that all available data were used for the calibration as suggested by Arsenault et al. [36]. The sensor at Necopastic is a GMON3 sensor, which is the same as the GMON instrument used by Choquette et al. [35], who estimated the measurement uncertainty to be in the order of 5-10% (for a SWE of less than 400 mm). In the Upper Yukon, a Campbell Scientific CS275 sensor was used at each station, with an uncertainty estimated to be in the order of ±15 mm (when the measured SWE is less than 300 mm), and ±15% otherwise.

Model Description: MASiN Snow Model
MASiN (Modèle Autonome de Simulation de la Neige) is a physically-based, multi-layer, snow model where the mass and energy balance for different snow layers are calculated while accounting for a maximum number of snow layers (MNSL) with a minimum snow layer depth of 1 cm (to conserve the stability of the iterative scheme). The MNSL value was set to 70 in the original version of the model in order to reduce the computational budget while conserving a certain inertia for energy and mass transfers in the snowpack. The energy balance takes into account the shortwave radiation according to the potential solar radiation theory of Lee [37] and is computed by considering the effect of cloud cover and vegetation as well as partitioning of the solar radiation into diffusive and direct shortwave radiations. The partitioning is done to model as accurately as possible the albedo and radiation absorption by the snowpack. The model also considers longwave radiation, sensible and latent turbulent heat fluxes, liquid water infiltration heat flux, and conduction flux. Table 3 summarizes the main physical processes simulated by MASiN as well as the related sub-processes and associated model parameters. Table 3. Physical processes simulated by MASiN and associated model parameters (N/A, stands for not applicable).

Shortwave radiation
Extraterrestrial irradiation N/A-

Effect of cloud and vegetation
Separation of direct and diffuse radiations k dir,min (Minimum ratio of direct shortwave radiation to total shortwave radiation) k dir,max (Maximum ratio of direct shortwave radiation to total shortwave radiation) Net shortwave radiation a dir,min (Minimum albedo for direct radiation) a diff,min (Minimum albedo for diffuse radiation)

Energy balances
Shortwave radiation β dir (Absorption coefficient for direct radiation) Due to the non-linearity of energy transfer between the snow layers, the internal computational time step was set to 30 s. While precipitation forcing is available at daily time steps, MASiN can only receive total precipitation values, or precipitation input partitioned into liquid and solid components. The case of the latter component, the model discretizes daily precipitation data into hourly time steps. For known rainfall and snowfall time series, rainfall data are discretized equally into 24 time steps, whereas snowfall amounts are discretized to as many hours as required while maintaining the minimum snow layer height (1 cm). If total daily precipitation data are only available, snowfall is assumed to occur when the hourly temperature is below 1 • C. To meet the minimum snow layer height requirement, snow can be redistributed into hourly time steps with temperature values below 1 • C. The temporal resolutions of the different input/output data are provided in Table 4. As previously mentioned, the objective of this study was to analyze the impact of the MNSL on snowpack (expressed in mm of SWE) modeling. The study was performed in two stages, as presented in Figure 2. In the first stage, MASiN was calibrated by specifying the MNSL. For each GMON station, a total of 21 calibrations were undertaken in order to simulate one to 20 MNSLs. The results were then compared against the full-extent configuration of MASiN (with a MNSL of 70). The rationale for this simulation was to reveal the existence of a threshold for the MNSL value, at which the modeling performance would not drop. This threshold was assumed to provide a limit for reducing the amount of interaction between the snow layers while avoiding any drop in SWE modeling performance. For this, the equifinality was studied to test the reliability of MASiN after decreasing the MNSL. This was achieved by selecting the top 10 best performances for each calibration of MASiN, at each GMON station, for all MNSL values.
In the second stage, the impact of the MNSL value on the MASiN's parameter set was analyzed based on the equifinality results obtained in the first stage. The equifinality analysis provides a range of values for different calibrated parameters for each value of the MNSL. Thus, the Pearson correlation coefficient between the 10 best calibrated parameters against the MNSL was calculated to estimate the possibility of a linear correlation between them. Here, the null hypothesis was defined as "there exists no correlation between each model parameter and the MNSL". When linear correlation exists, it would be possible to determine a more accurate range of variation for any future calibration where the MNSL value is modified. Table 5 presents the model parameters along with their upper and lower bounds and further highlights those that were kept for calibration as suggested in the literature [24]. To pair MASiN with HYDROTEL, daily total precipitation (Pt) values were partitioned into rainfall (R) and snowfall (S) using the following algorithm available in HYDROTEL, which is based on a temperature threshold (T thres ), the minimum (T min ), and maximum daily temperatures (T max ):

Calibration
It is noteworthy that the temperature threshold (T thres ) was calibrated while testing the performance of MASiN. As mentioned previously, in order to calibrate the MASiN parameters for a given MNSL, the dynamically dimensioned search (DDS) optimization algorithm was used. DDS can provide a set of calibrated parameter values, which is required to investigate the effect of equifinality. Based on the guidelines suggested by Tolson et al. [30], for each GMON station and each value of the MNSL, the calibration was grouped into 33 trials of 100 iterations and executed using a MATLAB DDS routine. Kling-Gupta efficiency (KGE) [38], given as follows, was then used as the objective function to quantify the goodness-of-fit between the measured and simulated SWEs: where µ s and µ o are the average values of the simulated (subscript s) and observed (subscript o) time series, respectively; σ s and σ o are the standard deviations of the simulated and observed time series, respectively; and r is the Pearson correlation coefficient. The larger is the KGE value, the more accurate the simulated series is, when compared to the observed series. At the end of each calibration, the top 10 best sets of parameter values were retained based on the final KGE value for each GMON station and the given MNSL value.

Results
In this section, the results of the study are presented in two separate subsections. Section 3.1 provides a comparison of the simulated SWE values as a function of MNSL. These results reveal whether any MNSL threshold value exists at which model performance could be maintained. Therefore, when comparing the results, the two major criteria of importance considered were SWE modeling accuracy and the required number of snow layers to simulate. The results are displayed using boxplots, where the MNSL values are shown on the abscissa, and the modeling performances are displayed on the ordinate. As the calibrated parameters directly affect different physical processes, studying their correlations with the MNSL values can provide valuable insights on how MASiN operates under different parameterizations. Thus, Section 3.2 investigates the correlations between the calibrated parameter values and MNSL values in order to identify how the latter can affect the individual physical processes modeled by MASiN when the MNSL decreases. This could ultimately provide information on how MASiN resolves the physical processes when the number of modeled snow layers varies.

Snow-Water Equivalent (SWE) Modeling
The top 10 best performances were compared against the MNSL values (see Figure 3). For each GMON station and all the simulation years pooled together, calibrations with a MNSL value of 1 provided a negative KGE value, and therefore are not shown in here.  The behavior of the model with respect to each station can also be examined. For Lower Fantail (Figure 3a), the full model configuration (i.e., with MNSL of 70) provided a median KGE value of 0.97. When reducing the MNSL value to 20 and down to nine, the modeling performance dropped only marginally from a median value of 0.92 to a median value of 0.89. When the MNSL was less than nine, the median KGE value varied between 0.87 and 0.89. Since all median KGE values were larger than 0.87, the best compromise to reduce the maximum number of simulated layers from a KGE value standpoint can be set to a value between two and nine for the GMON station at Lower Fantail. For the Lower Llewellyn station (Figure 3b), all median KGE values for any MNSL value including the fully-configured model (i.e., with 70 snow layers) were all in the range of 0.72 to 0.78. Therefore, it can be argued that from a KGE value point of view, a MNSL value of two can maintain the modeling performance at the same level for the full configuration of MASiN. For the Wheaton GMON station (Figure 3c), an almost similar profile to that of Lower Fantail was observed. The fully-configured model, however, provided a slightly higher median KGE of 0.97. The median KGE performances were achieved (0.92 to 0.96) for a range of nine to 20 layers. Below the MNSL value of nine, the median performance dropped to around 0.90 (MNSL values of five to seven), and then stabilized around 0.93 for the MNSL values of two to four. Therefore, any MNSL value between two and 20 can be considered to provide an acceptable modeling outcome. As shown in Figure 3d, for Necopastic, the fully-configured model was slightly outperformed by the model configuration with MNSL values between 15 and 20. Below this range of MNSL values, the modeling performance decreased slightly to reach a median performance of around 0.85 for a MNSL value of three. Similarly, for the Wheaton station, any MNSL value could be considered for a slight reduction in the modeling performance.
For the Lower Llewellyn, the relationship between KGE values and MNSL was quite different than those depicted at the other sites and counter intuitive from that viewpoint. In other words, while the relationships for the other sites could be characterized by a nearly parabolic shape, for Lower Llewellyn, it was almost linear with a negative slope. This suggests that there might be something not properly resolved with the GMON values.
To illustrate the uncertainties associated with each MNSL value,  Figure 3 were calculated using the whole calibration period of each GMON station, the KGE values provided in the captions of the aforementioned figures were solely calculated for the displayed winters. Thus, low KGE values, obtained for some parameter sets, can be explained by the compensation of the modeling for the other winters. Moreover, the negative KGE values for a MNSL of one can be explained by the underestimation of the modeled SWEs. In general, the larger the MNSL value, the better the simulation. This means that for the fully-configured model, the uncertainty is reduced, which is synonymous to an increasingly narrow equifinality. The timing of the dominant melting period was also improved, which was the best for the full model configuration (MNSL of 70).
As suggested earlier, the behavior of MASiN for the different configurations suggests that the observed SWEs at Lower Llewellyn might not be reliable, and overall not accurate enough to be considered for this study. However, they were kept here to demonstrate that a physically-based model can be quite useful to detect instrumentational errors or to infer other external factors that may corrupt observations.     Figure 8 shows, for each GMON station, the relative difference between the observed and simulated maximum SWE values for the top 10 best calibrated parameter sets for the calibration period for MNSL values of two, nine, and 70. Calibrations with MNSL values of nine (for Lower Fantail, Lower Llewellyn, and Necopastic) and two (for Lower Llewellyn) provided more accurate simulations of the maximum SWE. Thus, since the maximum SWE value reached over winter should lead to a more accurate annual hydrograph, these MNSL values can be referred to as thresholds. For the aforementioned GMON stations, and thus, the threshold values of the MNSL, the median values of the relative differences between the observed and simulated maximum SWE values were −4% at Lower Fantail, −13% at Lower Llewellyn, 2% at Wheaton, and 6% at Necopastic. Based on the above analysis for the GMON stations, a trade-off between the MNSL value and the SWE modeling accuracy could be obtained. In most cases, the fully-configured model with a MNSL value of 70 provided the best overall modeling performance as it could accurately outline the energy profile and mass transfer between the snow layers in the snowpack. Nevertheless, reducing the number of snow layers can still preserve a satisfactory level of SWE modeling performance in terms of the KGE. Depending on the required SWE modeling accuracy, any MNSL more than one can be considered as appropriate, even if a slight drop in performance can be observed at some GMON stations when the MNSL value is reduced below a certain threshold (e.g., MNSL of 9).

Influence of the Maximum Number of Snow Layer on the Calibrated Parameters
For each GMON station, the top 10 best performances provided a range of calibrated parameter values for each MNSL investigated. These sets of calibration parameter values were compared against the MNSL values in order to verify how they could influence the calibrated parameter values (refer to the Supplementary Document provided in the online version of this article). Thus, Pearson correlation coefficients between each calibrated parameter against the MNSL value were calculated and are shown in Table 6, along with the associated p-values. The p-values were compared against the significance level of 0.05. Table 6. Correlations of the top 10 best calibrated parameter sets (see Table 5 for the nomenclature) against the maximum number of snow layer (MNSL)value for each snow station when including all the snow stations (Global column). The values shown within the parentheses represent the p-values. Grey shaded cells indicate cases where the null hypothesis is not rejected, while blue shaded cells identify the parameters with significant Pearson correlation. Parameters with an absolute correlation superior to 0.3 were analyzed for their influence in MASiN.

Variable
Lower Globally, either the null hypothesis is not rejected for the correlations (shown as the shaded cases in Table 6), and so the linear correlation is not statistically different than zero, or the correlations are too low for describing a linear relationship between the parameters and the given MNSL. Consequently, the following discussion deliberated on the absolute correlations between the calibrated parameters and the MNSL values for correlations greater than 0.3.
The fresh snow minimum density (ρ ns ) showed a negative correlation for each GMON station when including all the stations in the linear regression. These correlations ranged from weak (Lower Llewellyn station) to moderate (Lower Fantail station), which indicates that the more snow layers are considered in MASiN, the less dense the fresh snow would be. This correlation is explained by the criteria for merging snow layers, which is based on the height of the layers [24]. As this criterion is met earlier in each winter for the smallest values of MNSL, MASiN considers that the height of the fresh snow layer is lower than the height of the fresh snow layer for the same meteorological conditions given the highest value of MNSL. Consequently, for the smallest values of MNSL, the fresh snow layer has a greater probability of being homogenized into older snow layers than staying distinct from the bottom layers.
For the settlement coefficient (K d ), there was a moderate positive correlation observed at Wheaton GMON station (R > 0.5). This positive correlation was also observed at Lower Fantail and Necopastic when including all the stations, but was weak (R > 0.2) at each of these stations. For Lower Llewellyn, the null hypothesis was not rejected. In MASiN, when the settlement coefficient increases, the height of the snow layer decreases. Thus, a positive correlation means that the compaction of the snow layers for the smallest values of MNSL is weaker than for the largest values of MNSL. So, by considering that the merging of the snow layer occurs earlier in each winter for the smallest values of MNSL, this compaction favors the homogenization of the snowpack from the top snow layers.
Finally, the maximum radiation coefficient (k SWmax ) had a weak negative correlation for the Lower Llewellyn GMON station (R < −0.3), while for other GMON stations and when including all the GMON stations, this negative correlation was weaker, but was not zero (R < −0.2). Increasing the maximum radiation coefficient in MASiN would increase the net direct shortwave radiation and the net diffuse shortwave radiation, which would favor the input of energy. This increase provides an additional amount of energy for each snow layer for the smallest values of MNSL compared to the largest ones, which can facilitate the melting process.
All in all, reducing the MNSL value would cause the fresh snow density and the maximum radiation coefficient to increase and the settlement coefficient to decrease. Accordingly, the melting process would be favored due to more input of energy from the shortwave radiation. Moreover, the top snow layers have a higher probability of being homogenized than the bottom layers. This is a phenomenon that is highlighted when the melting process induces a mass transfer between the snow layers, thus decreasing the layer heights when merging the snow layers with the smallest thickness. When compared to the acknowledged most sensitive phenomena reported in the literature (e.g., discrimination of rainfall and snowfall, albedo estimation, and the melting period), the correlation found for the radiation confirmed the significance of the melting period for an improved modeling performance. For MASiN, the influences of the fresh snow density and the settlement coefficient were not linked to a particular physical phenomenon, but are relevant for merging the smallest snow layers, and thus responsible for decreasing the number of snow layers to simulate, which was the motivation behind this study.

Discussion and Conclusions
Simulated SWE values are the main criteria when evaluating the performance of a snow model, whether it is used as a standalone module or as part of a hydrological model. Such an evaluation was undertaken in this study to ensure that a reasonable performance and modeling accuracy would be maintained by modifying the stratification of the modeled snowpack. In this study, the performance of the physically-based snow model MASiN was studied for a potential future integration in HYDROTEL. Modifying the MNSL was conceivable as it would decrease the modeling interactions within the snowpack. With this in mind, the objective of this study was to analyze the influence of the MNSL on model performance, namely SWE estimation.
The first part of the study compared the influence of the MNSL value on the SWE in terms of the KGE performance metric. The fully-configured MASiN (i.e., with a MNSL value of 70) provided the best overall performance, while reducing the MNSL to one caused the performance to drop significantly. However, globally speaking, using a MNSL value between two and 20 (rather than 70) would only marginally decrease the modeling performance to an acceptable level as the corresponding KGE values only dropped by less than 0.1 below that of the fully-configured model. In terms of the median value of the relative differences between the observed and simulated maximum SWEs, a slight performance drop could be observed for MNSL values less than or equal to nine. This was true for three GMON stations out of four analyzed in this study. Indeed, this study illustrated that a physically-based model can be quite useful to detect potential instrumental errors.
In the second part of the study, the impact of the MNSL value on the calibration parameter values was assessed. The analysis provided information on how the modeled physical processes behave when reducing the MNSL value. Thus, it becomes possible to adjust the MNSL value more adequately when the user assesses the need to amplify or reduce a modeled physical process to improve the SWE modeling. By taking into account each individual GMON station or by considering all of them together, it was shown that only some moderate correlation (|R| > 0.5) and weak correlation (|R| > 0.3) existed between the calibrated parameters and the MNSL values. Reducing the MNSL caused the fresh snow density and the maximum radiation coefficient to increase, while the settlement coefficient decreased. Consequently, by considering the influence of the MNSL on the fresh snow density and the settlement coefficient, MASiN preferentially reduces the height of the top snow layers. This means that it would preferentially homogenize the snow layers from the top of the snowpack since the snow layer merging condition is based on the minimum height, when compared to the fully configured version of the model. Moreover, the influence of a reduced MNSL value on the maximum radiation coefficient favors the melting process by adding more net shortwave radiation in the snowpack.
Finally, reducing the MNSL makes it possible to maintain a level of SWE modeling performance similar to that provided by the fully configured MASiN, when using the KGE performance metric. Indeed, although the modeling performances remained within an acceptable range, it was not possible to clearly identify the affected modeled physical processes. Meanwhile, the second part of this study showed that reducing the MNSL did affect a few model parameters, allowing the identification of the modeled processes influenced by the change in the MNSL values. Consequently, the losses in modeling accuracy were primarily associated with snow inputs (i.e., fresh snow density), snow layer settlement, and melting process (i.e., maximum amount of radiation). While snow inputs and melting process have already been identified as significant processes in previous studies ( [18,[20][21][22]), this study showed that the settlement process was identified as an additional phenomenon affected by the reduction of the MNSL values. However, the influence on the settlement process was related in all likelihood to the snow layer merging conditions, which are specific to MASiN. Moreover, Domine et al. [16,17] noticed that the modeling accuracy of the vertical profile of snow density, which are associated with the snowpack thermal properties, was paramount in the modeling of the groundwater budget. Indeed, the vertical profile of snow density simulated by MASiN was directly affected by the fresh snow density and the settlement coefficient; the latter parameters were influenced by a reduction of the MNSL values, causing a drop in modeling performance. Meanwhile, before pairing MASiN with HYDROTEL, first, it is important to consider a methodology to spatially extrapolate some of the input data like hourly relative humidity and wind speed. Once a methodology is developed to spatially extrapolate the input data for MASiN, and after pairing the model with HYDROTEL, the runoff modeling accuracy can be estimated. Finally, the framework introduced in this paper has the potential to be applied to other physically-based snow models that provide a means to adjust the number of simulated snow layers and as long as possible to save the model results for comparison purposes.
Supplementary Materials: The following information can be found online at http://www.mdpi.com/2073-4441/ 12/12/3449/s1, Figure S1: Snow layer density triggering the metamorphism phenomenon of the snow layer against the MNSL at GMON LF station, Figure S2: Snow layer density triggering the metamorphism phenomenon of the snow layer against the MNSL at GMON LL station, Figure S3: Snow layer density triggering the metamorphism phenomenon of the snow layer against the MNSL at GMON W station, Figure S4: Snow layer density triggering the metamorphism phenomenon of the snow layer against the MNSL at GMON Neco station, Figure S5: Fresh snow minimum density against the MNSL at GMON LF station, Figure S6: Fresh snow minimum density against the MNSL at GMON LL station, Figure S7: Fresh snow minimum density against the MNSL at GMON W station, Figure S8: Fresh snow minimum density against the MNSL at GMON Neco station, Figure S9: Maximum retention capacity of the snow layer against the MNSL at GMON LF station, Figure S10: Maximum retention capacity of the snow layer against the MNSL at GMON LL station, Figure S11: Maximum retention capacity of the snow layer against the MNSL at GMON W station, Figure S12: Maximum retention capacity of the snow layer against the MNSL at GMON Neco station, Figure S13: Settlement coefficient against the MNSL at GMON LF station, Figure S14: Settlement coefficient against the MNSL at GMON LL station, Figure S15: Settlement coefficient against the MNSL at GMON W station, Figure S16: Settlement coefficient against the MNSL at GMON Neco station, Figure S17: Ground heat flux against the MNSL at GMON LF station, Figure S18: Ground heat flux against the MNSL at GMON LL station, Figure S19: Ground heat flux against the MNSL at GMON W station, Figure S20: Ground heat flux against the MNSL at GMON Neco station, Figure S21: Atmospheric temperature threshold associated with the fresh snow minimum density against the MNSL at GMON LF station, Figure S22: Atmospheric temperature threshold associated with the fresh snow minimum density against the MNSL at GMON LL station, Figure S23: Atmospheric temperature threshold associated with the fresh snow minimum density against the MNSL at GMON W station, Figure S24: Atmospheric temperature threshold associated with the fresh snow minimum density against the MNSL at GMON Neco station, Figure S25: Snow cover surface roughness against the MNSL at GMON LF station, Figure S26: Snow cover surface roughness against the MNSL at GMON LL station, Figure S27: Snow cover surface roughness against the MNSL at GMON W station, Figure S28: Snow cover surface roughness against the MNSL at GMON Neco station, Figure S29: Reduction coefficient of the turbulent trade against the MNSL at GMON LF station, Figure S30: Reduction coefficient of the turbulent trade against the MNSL at GMON LL station, Figure S31: Reduction coefficient of the turbulent trade against the MNSL at GMON W station, Figure S32: Reduction coefficient of the turbulent trade against the MNSL at GMON Neco station, Figure S33: Minimum radiation coefficient against the MNSL at GMON LF station, Figure S34: Minimum radiation coefficient against the MNSL at GMON LL station, Figure S35. Minimum radiation coefficient against the MNSL at GMON W station, Figure S36: Minimum radiation coefficient against the MNSL at GMON Neco station, Figure S37: Maximum radiation coefficient against the MNSL at GMON LF station, Figure S38: Maximum radiation coefficient against the MNSL at GMON LL station, Figure S39: Maximum radiation coefficient against the MNSL at GMON W station, Figure S40: Maximum radiation coefficient against the MNSL at GMON Neco station, Figure S41: Threshold temperature of precipitation separation against the MNSL at GMON LF station, Figure S42: Threshold temperature of precipitation separation against the MNSL at GMON LL station, Figure S43: Threshold temperature of precipitation separation against the MNSL at GMON W station, Figure S44: Threshold temperature of precipitation separation against the MNSL at GMON Neco station.