Seasonal Variability of Snow Density in the Spanish Pyrenees

: Spanish latitudes and meteorological conditions cause the snow phenomena to mainly take place in mountainous areas, playing a key role in water resource management, with the Pyrenees as one of the most important and best monitored areas. Based on the most signiﬁcant dataset of snow density (SDEN) in the Spanish Pyrenees for on-site manual samples and automatic measurements, in this study, single and multiple linear regression models are evaluated that relate SDEN with intra-annual time dependence and other drivers such as the seasonal accumulated precipitation, 7-day average temperatures, snow depth (SD) and elevation. The seasonal accumulated precipitation presented a more dominant inﬂuence than daily precipitation, usually being the second most dominant SDEN driver, followed by temperature. Average temperatures showed the best ﬁtting to SDEN. The results showed similar densiﬁcation rates ranging widely from 0.7 × 10 3 kg/L/day to 2 × 10 3 kg/L/day without showing a spatial pattern. The densiﬁcation rate for the set of manual samples was set to 1.2 kg/L/day, very similar to the set of automatic measurements (1.3 kg/L/day). The results increase knowledge on SDEN in the Pyrenees. The SDEN regression models that are given in this work may allow us, in the future, to estimate SDEN, and consequently Snow Water Equivalent (SWE), using an economical and extensive SD and meteorological network, although the high spatial variability that has been found must be regarded. Estimating a relationship between SDEN and several climate drivers enables us to take into account the impact of climate variability on SDEN.


Introduction
Snow in Spain plays a key role in water resource management and occurs essentially in mountainous areas due to Spanish latitudes and meteorological conditions [1][2][3][4][5]. One of the most important Spanish mountain regions is the Pyrenees, which extends about 490 km between Spain and France, connecting the Atlantic Ocean at its eastern limit with the Mediterranean Sea at its western limit (both of which influence the climate of the Pyrenees), and with Aneto Peak as its highest point (3404 m altitude). From the different parameters defining snow-related phenomena, snow density (SDEN) relates snow depths (SDs) with the snow water equivalent (SWE). Since SD measurements are much more frequent than the available SWE data, better knowledge and parameterization of SDEN will enable the achievement of a more accurate estimation of the SWE based on SD measurements. Moreover, as wind data become more available, hydrological models can determine a more accurate estimation of forced-convection snowmelt and snow transport, for which SDEN is needed. However, SDEN is a complex parameter that can vary temporally, spatially and even within the snowpack profile in the vertical direction.
The density of freshly fallen snow depends on three main climatic conditions: temperature, wind and humidity. Typical proposed density values for recently fallen snow are found between 0.07 and 0.15 kg/L, with 0.1 kg/L being the most common value [6][7][8].
the SWE estimated from the computed SDEN was less than 15% and was similar to that obtained by relating SDEN directly to SD [40].
The Spanish Pyrenees hold more than one hundred snow poles and up to thirteen automatic snow monitoring devices, or telenivometers (TNMs), operating since 2008. This instrumentation allows for the study of the seasonal variability of snowfall events [41]. Additionally, the Spanish and Ebro Water Authority have been carrying out on-site SDEN measurements since the early 1990s through the ERHIN program ("Estudio de Recursos Hídricos procedentes de la INnivación"-Study of Snow Water Resources) [42]. TNM data show significant variability in annual SD. Western regions (Atlantic Ocean weather) give the most annual accumulation in the winter months, while eastern regions (Mediterranean weather types) show more homogenous accumulation over winter, spring and late autumn. It is interesting to note that western TNMs also show behaviours similar to that observed for the western coast of the US, which is characterized by plenty of precipitation and relatively high air temperatures, both possible contributors to early season high SDEN [43].
The novelty of this paper remains in the definition of SDEN variability in the Spanish Pyrenees, using the vast amount of data from the ERHIN program obtained from manual sampling and non-destructive automatic measurements. A series of relationships are established with other parameters, and its temporal and spatial variability is characterized. Few comparable works in the world have used such a wide database (SNOTEL or WSL Institute for Snow and Avalanche Research SLF), and there are no similar precedents in the study area. To achieve this goal, a series of single and multiple linear regressions are conducted, relating SDEN with the time evolution variability, seasonal accumulated precipitation, SD, elevation and average temperatures.

Materials and Methods
The ERHIN program [44] started in 1986, studying snow phenomena in river basins in the mountainous regions of the Pyrenees, Sierra Nevada, the Central System and the Cantabrian Mountains. In this framework, on-site manual SDEN measurements were carried out and a dense network of snow poles and TNMs (SWE and SD non-destructive measurements) were installed. The Pyrenees is the densest region monitored, with 13 TNMs and more than 100 snow poles ( Figure 1). The TNMs cover a range about 200 km long in the Pyrenees and are numbered from west to east.

On-Site Manual Sampling and Measurements
Since 1987, on-site manual SDEN measurements have been carried out at certain snow poles. At the same time, SD measurements were taken for the whole snow pole network. Sampling elevations range from 1440 m to 2615 m and cover up to three seasonal campaigns between December and May. Figure 2 shows a 0.8 m stackable sampling tube of constant 53.8 mm diameter inside a section (22.72 cm 2 ) and 2.5 mm thick, with up to 4 m of maximum penetration depth. This SWE tube was handmade in 1986 at the Technical University of Valencia. Knowing the weight and dimensions of the core, and therefore its volume, by directly weighing it on a precision weighing scale, the net weight of the snow core is obtained, and the density is calculated.

On-Site Manual Sampling and Measurements
Since 1987, on-site manual SDEN measurements have been carried out at certain snow poles. At the same time, SD measurements were taken for the whole snow pole network. Sampling elevations range from 1440 m to 2615 m and cover up to three seasonal campaigns between December and May. Figure 2 shows a 0.8 m stackable sampling tube of constant 53.8 mm diameter inside a section (22.72 cm 2 ) and 2.5 mm thick, with up to 4 m of maximum penetration depth. This SWE tube was handmade in 1986 at the Technical University of Valencia. Knowing the weight and dimensions of the core, and therefore its volume, by directly weighing it on a precision weighing scale, the net weight of the snow core is obtained, and the density is calculated.  Table 1 shows the on-site manual sampling and elevation for sites in the same locations as 11 TNM sites. The number of samples is also included. The complete list of 35 snow pole sites and the number of samples are given as supplementary data. Figure 3 shows the elevation distribution for all of the samples taken, ranging from 1810 m to 2615 m for snow pole sites.   Table 1 shows the on-site manual sampling and elevation for sites in the same locations as 11 TNM sites. The number of samples is also included. The complete list of 35 snow pole sites and the number of samples are given as supplementary data. Figure 3 shows the elevation distribution for all of the samples taken, ranging from 1810 m to 2615 m for snow pole sites.   (1) Table 2 shows the name, situation, elevation, number of daily samples and sampling period for each TNM, which were located above 1800 m elevation and with 2600 m as the maximum elevation. It should be noted that SDEN calculated according to Equation (1) is not sensitive to a higher time resolution. Figure 1 shows the situation of the TNMs in the Pyrenees, (named from west to east), while Figure 5 shows the elevation distribution for the automatic TNM measurements. It is interesting to mention that the Eriste TNM is located at approximately the same distance from the Atlantic Ocean and the Mediterranean Sea.

On-Site Non-Destructive Measurements
Although the TNMs have been operating for a shorter period than the period covered by manual samples, daily data were registered, with more than 20,000 samples and officially provided by the Ebro Water Authority (http://www.saihebro.com/saihebro/index.php?url=/datos/usos/mapa:TNG/tipoestacion:TN accessed on 4 June 2021). Additionally, the maximum, average and minimum daily temperatures were available.

On-Site Non-Destructive Measurements
First Cosmic-Ray Neutron (CRN) attenuation telenivometers (TNMs) in 2008 along the Spanish Pyrenees (Figure 4), and two more were recen 2015 and 2018. TNMs allow the obtaining of SWE high time resolution dat SD vertical data (1 cm resolution). From those values, SDEN can be estim to Equation (1): = Table 2 shows the name, situation, elevation, number of daily samples period for each TNM, which were located above 1800 m elevation and with maximum elevation. It should be noted that SDEN calculated according to not sensitive to a higher time resolution. Figure 1 shows the situation of th Pyrenees, (named from west to east), while Figure 5 shows the elevation d the automatic TNM measurements. It is interesting to mention that the Eri cated at approximately the same distance from the Atlantic Ocean and the Sea.
Although the TNMs have been operating for a shorter period than the p by manual samples, daily data were registered, with more than 20,000 sam cially provided by the Ebro Water Authority (http://www.saihebro.co dex.php?url=/datos/usos/mapa:TNG/tipoestacion:TN accessed on 4 June 20 ally, the maximum, average and minimum daily temperatures were availa  Table 2 shows the name, situation, elevation, number of daily samples and sampling period for each TNM, which were located above 1800 m elevation and with 2600 m as the maximum elevation. It should be noted that SDEN calculated according to Equation (1) is not sensitive to a higher time resolution. Figure 1 shows the situation of the TNMs in the Pyrenees, (named from west to east), while Figure 5 shows the elevation distribution for the automatic TNM measurements. It is interesting to mention that the Eriste TNM is located at approximately the same distance from the Atlantic Ocean and the Mediterranean Sea.

Other Data and Sources of Error
Supplementary Materials are provided for both on-site manual sampling and on-site non-destructive (TNM) daily measurements.
The daily data for precipitation from the closest meteorological stations within the same period as the TNM measurements are also given in the Supplementary Material.
The sources of error and uncertainties associated with snow sampling include the loss of part of the sample due to snowpack collapse when the sampler encounters hard layers. When the sampler is extracted from the snowpack for weighing, snow can also be lost [45] or stuck to the outside of the sampler. These errors may be increased by a snowpack with ice layers that give the false perception of reaching the ground, basal ice and layers of depth hoar and non-cohesive crystals. New snow or wet snow can also increase sampling errors, especially if the observer is collecting core samples to be weighed later rather than using a spring balance and tube cradle. The experience of the observer plays Although the TNMs have been operating for a shorter period than the period covered by manual samples, daily data were registered, with more than 20,000 samples and officially provided by the Ebro Water Authority (http://www.saihebro.com/saihebro/index.php? url=/datos/usos/mapa:TNG/tipoestacion:TN accessed on 4 June 2021). Additionally, the maximum, average and minimum daily temperatures were available.

Other Data and Sources of Error
Supplementary Materials are provided for both on-site manual sampling and on-site non-destructive (TNM) daily measurements.
The daily data for precipitation from the closest meteorological stations within the same period as the TNM measurements are also given in the Supplementary Material.
The sources of error and uncertainties associated with snow sampling include the loss of part of the sample due to snowpack collapse when the sampler encounters hard layers. When the sampler is extracted from the snowpack for weighing, snow can also be lost [45] or stuck to the outside of the sampler. These errors may be increased by a snowpack with ice layers that give the false perception of reaching the ground, basal ice and layers of depth hoar and non-cohesive crystals. New snow or wet snow can also increase sampling errors, especially if the observer is collecting core samples to be weighed later rather than using a spring balance and tube cradle. The experience of the observer plays an important role in reducing potential errors due to human failure. Depending on the temperature, new and wet snow will tend to stick in the tube, and this will result in an underestimation of the SWE, which could exceed 10%. Smaller diameter cutters (down to 20 cm 2 ) were more prone to plugging as they encountered ice lenses, and they were more likely to induce the collapse of non-cohesive layers under the cutter, resulting in an underestimation of the total SWE [46]. The ideal cutter area was about 30 cm 2 , demonstrated by the low error percentage of the ESC30 sampler, which ranged from a 5% overestimation to a 2% underestimation [47]. The uncertainty of density measurements in non-ideal snow conditions is approximately within 10 to 15% [48]. As a control measure, random validation tests were carried out for manual samples and measurement protocols were followed [46].
Errors in TNM measurements can occur due to non-environmental issues such as instrument malfunction and incorrect instrument calibration (or calibration drift), CRN uncertainties due to changes in mid-season soil moisture levels and the undervaluation of precipitation and gaps or false null data from other meteorological variables. In this regard, TNM data are periodically calibrated with manual sampling during late snow conditions [42].

Analysis Methodology
A previous screening process was carried out, with the identification of outliers (a range check for reasonable values between zero and the maximum possible SDEN for the site), gaps and null data, discarding sites with little information or information that was not representative.
TNM measurements are fully automated and, even though they are frequently calibrated, are not quality controlled. Thus, the calculated SDEN values that significantly exceeded the expected maximum and minimum values at each site were removed. The number of samples for each TNM is given after quality-control procedures. All negative and null values were eliminated from the records, for both the SWE and SD, which are necessary to estimate SDEN. Consequently, TNMs numbers 8, 14 and 15 (Salenques, Sarrios-Formigal and Besurta) were discarded.
Manual samples, much less prone to uncertainties [37], underwent a similar screening process. On-site manual sample sites with enough data in the same location as TNM sites were selected and pooled. The previous screening process reduced 35 on-site manual sites with 841 samples to five on-site manual sites matching TNM sites with 284 samples. SDEN statistics were calculated for both manual samples and automatic measurements during the entire year (annual period), as well as for the accumulation period (winter snow, until 21 March) and for the melting period (spring snow, from 21 March).
Automatic TNM sites with a significant correlation between SDEN and the day of the year (time evolution) were selected and grouped (the set of TNMs) for automatic TNM measurements, where day 1 corresponds to 1 October. The results were compared with the SDEN-day of the year correlation for manual samples.
Independent variables were analysed for the selected automatic TNM sites. The dominant climatological drivers of annual SDEN variability were extracted by establishing single linear regressions between each of the predictor terms and SDEN. The selected climate and spatial predictors include the time evolution, precipitation, accumulated precipitation, SD and average temperature.
Multiple linear regressions (MLRs) between SDEN and these predictors were studied for automatic measurements and compared with manual samples. The effect of elevation on the grouped series (the set of TNMs) was evaluated. Predictors were added progressively in the MLR model from more to less dominant with a significance threshold of 0.05.

Results
In this section, the SDEN results are analysed, identifying dominant climatological variables and studying MLR models for SDEN. Results from both the manual and automatic measurements are presented. Table 3 shows a summary of the SDEN statistics obtained by on-site manual data for two annual periods, winter or early season (until 21 March) and spring or late season (21 March-30 June). The word "set" refers to all grouped data. The average SDEN and coefficient of variation (C.V.) were very similar among sites for every period, without spatial patterns from east to west locations. Renclusa showed the lowest average SDEN and the highest C.V. due to a very low winter average SDEN and very marked differences between early and late season snow. Bachimaña and Ordiceto showed the highest average SDEN. As expected, the average SDEN value was greater for the late season (spring) than for the early season (winter). The C.V. is small for all sites studied and the set of samples, being the greatest during the early season (winter). Table 4 shows a summary of the SDEN statistics obtained from TNM data for two annual periods, winter or early season (until 21 March) and spring or late season (21 March-30 June). The set of TNMs refers to TNMs numbers 1, 2, 4, 6, 7 and 9. As for manual samples, except for Eriste, the TNM data showed very similar SDEN averages among sites for every period, without spatial patterns from east to west locations. Eriste showed the highest average SDEN, while Canal Roya's average SDEN was the lowest; the C.V. is again small, with all TNMs below 5%, but there were some differences among sites, since the amount of data between them is variable. The Ordiceto TNM showed the greatest C.V. Once again, the C.V. was usually greater for the early season (winter) than for the late season (spring).

SDEN Statistics
A comparison of statistics between manual and automatic data (see Table 5) showed small differences in terms of the average SDEN, although the C.V. values are higher for the automatic data. The set of manual samples and the set of TNM data show small differences for the average SDEN, especially for the annual and early-season periods. Renclusa showed the lowest annual SDEN average for both the manual and automatic data, while the highest annual average SDEN was represented by different locations for each technique.

Most Representative Automatic TNM Sites
The most representative automatic TNM sites were selected and grouped (the set of TNMs) for the automatic TNM data, with a significant correlation between SDEN and the day of the year (intra-annual time dependence), where day 1 corresponds to 1 October. An intra-annual time dependence for SDEN was found using Equation (2): The coefficients A and B are given in Table 6 along with the site elevation of each TNM and the average elevation of the set of TNMs, with the coefficient of determination (R 2 ) being computed in each case. The relationship found is important as it describes SDEN when data from other climatic variables are not available. The coefficient A represents the densification rate, ranging from 0.7 × 10 −3 kg/L/day to 2 × 10 −3 kg/L/day; the coefficient B represents the initial SDEN, ranging from 0.1 kg/L to 0.3 kg/L. Neither of the two parameters shows a spatial pattern. Automatic TNM sites with a significant coefficient of determination (R 2 ) were selected and grouped (the set of TNMs) for automatic TNM measurements. The densification rate for the data from the set of TNMs is set to 1.3 × 10 −3 kg/L/day, with an initial SDEN of 0.196 kg/L. Ordiceto, placed in the central Pyrenees, showed the highest densification rate and coefficient of determination (R 2 ) and the lowest initial SDEN. Aixeus, in the eastern part of the Pyrenees, showed the lowest densification rate and coefficient of determination (R 2 ), while the highest initial SDEN was found in Eriste. Canal Roya, Bachimaña, Lapazosa, Eriste, Airoto and Aixeus showed low densification rates and high initial SDENs, while Quimboa, Izas, Ordiceto and Renclusa showed the opposite behaviour. Therefore, no spatial pattern can be observed.
Regarding (R 2 ), TNMs numbers 1, 2, 4, 6, 7 and 9 are selected. TNM number 9 (Eriste) is selected as the most representative of the eastern part of the Pyrenees. Figure 6 shows the SDEN intra-annual time dependence for the selected TNMs and the trend line for the set of TNMs.   Table 7 shows the values of the coefficients A and B for manual sampling in the most representative TNM sites and the set of manual samples, which enables the establishment of the intra-annual time dependence of SDEN using Equation (2). The average elevation for the set of manual samples and the coefficient of determination (R 2 ) are shown, too. For manual samples, the densification rates and an initial SDEN range have less variability than the automatic measurements, from 1.1·10 −3 kg/l/day to 1.8·10 −3 kg/l/day and from 0.13 kg/l to 0.19 kg/l, respectively. The densification rate from the set of manual samples is set to 1.2·10 −3 kg/l/day, with an initial SDEN of 0.189 kg/l, showing values very similar to but lower than the values of the automatic measurements. Common sites for both manual and TNM data are Bachimaña, Ordiceto and Renclusa. Bachimaña showed the same coefficient of determination but a higher densification rate and a lower initial SDEN for manual samples than for automatic data. Ordiceto showed a lower coefficient of determination and densification rate, but a higher initial SDEN for the manual samples than for the automatic data. Renclusa had a higher coefficient of determination and the same densification rate, but a lower initial SDEN.  Table 7 shows the values of the coefficients A and B for manual sampling in the most representative TNM sites and the set of manual samples, which enables the establishment of the intra-annual time dependence of SDEN using Equation (2). The average elevation for the set of manual samples and the coefficient of determination (R 2 ) are shown, too. For manual samples, the densification rates and an initial SDEN range have less variability than the automatic measurements, from 1.1 × 10 −3 kg/L/day to 1.8 × 10 −3 kg/L/day and from 0.13 kg/L to 0.19 kg/L, respectively. The densification rate from the set of manual samples is set to 1.2 × 10 −3 kg/L/day, with an initial SDEN of 0.189 kg/L, showing values very similar to but lower than the values of the automatic measurements. Common sites for both manual and TNM data are Bachimaña, Ordiceto and Renclusa. Bachimaña showed the same coefficient of determination but a higher densification rate and a lower initial SDEN for manual samples than for automatic data. Ordiceto showed a lower coefficient of determination and densification rate, but a higher initial SDEN for the manual samples than for the automatic data. Renclusa had a higher coefficient of determination and the same densification rate, but a lower initial SDEN.

Identification of Dominant Variables
In addition to the time dependence, the selected predictors included the SD (cm), seasonal accumulated daily precipitation (PPacum (mm)) and 7-day average temperature (Tave7d) (°C/7 days). The seasonal accumulated precipitation was taken into account from the first day that snow accumulation starts and presented a more dominant influence than the daily precipitation, as SDEN depends on historical intra-annual evolution. The 7-day average temperatures have shown the best fitting to SDEN. SDEN depends on punctual temperatures for the accumulation process, but during melting cycles, average temperatures may describe it better [33]. The effects of the maximum and average daily temperatures were similar, while the minimum daily temperature was less dominant. Average daily temperatures were chosen, as they may be more available. Table 8 shows the correlation values, with the intra-annual time dependence and seasonal accumulated precipitation being the most dominant variables.

Identification of Dominant Variables
In addition to the time dependence, the selected predictors included the SD (cm), seasonal accumulated daily precipitation (PPacum (mm)) and 7-day average temperature (Tave7d) ( • C/7 days). The seasonal accumulated precipitation was taken into account from the first day that snow accumulation starts and presented a more dominant influence than the daily precipitation, as SDEN depends on historical intra-annual evolution. The 7-day average temperatures have shown the best fitting to SDEN. SDEN depends on punctual temperatures for the accumulation process, but during melting cycles, average temperatures may describe it better [33]. The effects of the maximum and average daily temperatures were similar, while the minimum daily temperature was less dominant. Average daily temperatures were chosen, as they may be more available. Table 8 shows the correlation values, with the intra-annual time dependence and seasonal accumulated precipitation being the most dominant variables. For the Quimboa TNM, the PPacum driver was the most dominant, while Tave7d was the most dominant for the Renclusa TNM. For the rest of the selected TNMs (Izas, Bachimaña, Ordiceto and Eriste), the intra-annual time dependence was the main variable.
The SDEN linear correlation with SD is almost negligible, except for the Ordiceto and Eriste TNMs. However, it may give useful information, as the significance threshold of SDEN is reduced for high SD values. Figure 8 represents SDEN's variability on SD for the set of TNMs, divided into an early-mid season snowpack and mid-late season snowpack. As observed, SDEN has a wide range of values with low SD, especially for the early season snowpack, reducing SDEN variability as SD increases. For the Quimboa TNM, the PPacum driver was the most dominant, while Tave7d was the most dominant for the Renclusa TNM. For the rest of the selected TNMs (Izas, Bachimaña, Ordiceto and Eriste), the intra-annual time dependence was the main variable.
The SDEN linear correlation with SD is almost negligible, except for the Ordiceto and Eriste TNMs. However, it may give useful information, as the significance threshold of SDEN is reduced for high SD values. Figure 8 represents SDEN's variability on SD for the set of TNMs, divided into an early-mid season snowpack and mid-late season snowpack. As observed, SDEN has a wide range of values with low SD, especially for the early season snowpack, reducing SDEN variability as SD increases.

Multiple Linear Regressions (MLRs)
The MLRs between SDEN and the evaluated predictors were studied for the automatic measurements and compared with the manual samples. Those MLRs follow Equation (3): The coefficients Ai are given in Table 9 for each site studied. As mentioned before, predictors were added progressively in the MLR model from more to less dominant, indicating the incrementally adjusted correlation coefficient (adjusted R 2 ). Moreover, the effect of elevation for the data from the set of TNMs was also evaluated as a driver, thus the MLR followed Equation (4)

Multiple Linear Regressions (MLRs)
The MLRs between SDEN and the evaluated predictors were studied for the automatic measurements and compared with the manual samples. Those MLRs follow Equation (3): The coefficients A i are given in Table 9 for each site studied. As mentioned before, predictors were added progressively in the MLR model from more to less dominant, indicating the incrementally adjusted correlation coefficient (adjusted R 2 ). Moreover, the effect of elevation for the data from the set of TNMs was also evaluated as a driver, thus the MLR followed Equation (4):  (3) and (4).

Name
A 1 (day) Therefore, (R 2 ) for A 4 (SD) represents the adjusted correlation coefficient for the MLR model and Equation (3), while (R 2 ) for A 5 (Elevation) represents the adjusted correlation coefficient for the MLR model and Equation (4).
The MLR models improve the correlation compared with the single linear regression for intra-annual temporal dependence (A 1 ), especially for those with a low adjusted R 2 , such as the Renclusa and Eriste TNMs. The adjusted R 2 ranges from 0.54 for the Bachimaña TNM to 0.83 for the Ordiceto TNM. The seasonal accumulated precipitation (PPacum) improves the correlation for the Quimboa and Ordiceto TNMs. The 7-day average temperature (Tave7d) plays a key role in describing SDEN for the Renclusa TNM and influences SDEN moderately for the Izas and Eriste TNMs. SD had an important effect in describing SDEN for the Bachimaña, Renclusa and especially the Eriste TNMs, and consequently, in the overall set of TNMs. The set of TNMs shows an adjusted R 2 of 0.57. However, adding elevation showed very little improvement, as the SD and elevation had a great collinearity.

Discussion and Conclusions
The temporal and spatial variability of snow density (SDEN) were studied in this paper, using the biggest SDEN data bank for on-site manual samples and automatic measurements in the Spanish Pyrenees. More than 375 manual samples and 21,000 automatic TNM data from the ERHIN program and Ebro Water Authority were used to define SDEN variability in the Spanish Pyrenees. Single and multiple linear regressions were conducted, relating SDEN with the time dependence and other drivers such as the seasonal accumulated precipitation, average temperatures, snow depth (SD) and elevation.
The automatic measurements (TNMs) provided a better description of SDEN variability than the manual SDEN sampling, as the data volume was about 20 times larger.
For both the manual samples and the TNM data, the average SDEN values and the C.V. were very similar, not following a spatial pattern. This similarity was proven as well for the set of manual samples and the set of TNM measurements. Additionally, the C.V. for the winter snow was usually larger than that of the spring snow. The Izas TNM showed similar SDEN values compared with previous works, with 0.2 and 0.3 kg/L during the first months and increasing to 0.6 kg/L at the end of the season [30]. However, it should be noted that similar values for the average and the C.V. are not representative of SDEN variability, either temporally or spatially. Intra-annual time dependence was shown to be the predominant driver for almost all TNM sites. The densification rates ranged widely from 0.7 × 10 −3 kg/L/day to 2 × 10 −3 kg/L/day, without showing a spatial pattern, being 1.3 × 10 −3 kg/L/day for the data from the set of TMNs. The densification rate for the set of manual samples was set to 1.2 kg/L/day, which is very similar to the automatic measurements. The densification rates were higher than those estimated for alpine regions in the former Soviet Union and the US, with average spring SDENs (kg/L) and snow densification rates (kg/L/day) of 0.25 kg/L and 0.79 × 10 −3 kg/L/day, respectively, for the former Soviet Union and 0.31 kg/L and 1.07 × 10 −3 kg/L/day, respectively, for the US. [37].
The seasonal accumulated precipitation presented a more dominant influence than daily precipitation, being the second most dominant SDEN driver for all sites except two. Historical behaviour showed a greater influence than daily behaviour. The third most important driver was the temperature. The 7-day average temperatures showed the best fitting to SDEN. Snow densification depends on punctual temperatures for the accumulation process, but during melting cycles, average temperatures may describe it better [33]. SD was found to be less significant. For low SD, both low and high SDEN can be found, depending on an early season dominated by accumulation processes or a late season dominated by melting processes [33,36]. For the early season, SDEN will depend on the proportion of precipitation falling as rain and snow, so it is also possible to reach high SDEN values. However, the significance threshold of SDEN is distinctly reduced for high SD values. SDEN for the late-season snow shows higher values than the early-season snow for maximum SD as densification continues, even though the SD may keep similar values for a certain period.
The established significance order for the drivers matched with previous studies, where precipitation was the dominant climate variable at most sites, followed by the average temperature and melt-refreeze (MRF) events [37].
Multiple Linear Regression (MLR) models were created for the automatic TNM measurements and were compared with the manual samples. The predictors were added progressively from more to less dominant. Elevation was added as a predictor for the set of TMNs series, showing that the SD was enough to describe that component of SDEN variability. The adjusted R 2 ranged from 0.54 for the Bachimaña TNM to 0.83 for the Ordiceto TNM, and the set of TNMs showed an adjusted R 2 of 0.57. These correlation coefficients were higher than the adjusted R 2 found for single linear regressions based on its time evolution, both for the manual samples and the automatic TNM data, with the adjusted R 2 ranging from 0.27 for Eriste to 0.73 for Ordiceto and 0.44 for the set of TNMs. This greater correlation means the existence of variability in SDEN can be explained by meteorological drivers, which are highly variable in time and space.
The R 2 values are similar to those obtained in the central Spanish Pyrenees, although the scale of work was different [40], and to those for spring snow in Alpine regions of the former Soviet Union and the US (R 2 = 0.68) [37]. These R 2 were obtained using MLR models with five variables (Precipitation, interactive SDmax * Temp, Cooling Degree Day, Latitude and Elevation). In the Swiss Alps, according to three elevation range classes (≥2000 m, ≥1400 and <2000 m, <1400 m) and 12 seasonal classes (months), SDEN was fitted to the SD with the R 2 values ranging from 0.02 to a maximum of 0.58 [36].
The results presented in Table 6 show that, in order to increase the available data and reduce uncertainties to enlarge the knowledge of SDEN in the Spanish Pyrenees, further revision and calibration should be carried out on TNMs 3 (Canal Roya), 5 (Lapazosa), 10 (Airoto) and 11 (Aixeus). The last two are essential to represent spatial variability due to the Mediterranean climate. Other TNMs, such as numbers 8 (Salenques), 14 (Sarrios-Formigal) and 15 (Besurta), should enlarge their data bank to be representative and contribute to similar studies. Good maintenance of the meteorological network is crucial for reducing outliers, gaps or false nulls and the undervaluation of precipitation.
The validation of these SDEN regressions could be done in other Spanish mountainous areas or other regions with similar or different climates. Results regarding the eastern locations should fit better for the Atlantic Climate (the Quimboa TNM), while those regarding western locations should fit better for the Mediterranean Climate (Eriste), although no spatial patterns were observed.
The results contribute to enlarging the knowledge of SDEN in the Pyrenees and enable it to be described in a more accurate way. The SDEN regression models that are given in this work may allow us in the future to estimate SDEN, and consequently SWE, using an economical and extensive SD and meteorological network, although the high spatial variability that has been found must be taken into account. Additionally, these SDEN regression models can be implemented in hydrological snow models to describe more complex snow transport or forced-convection phenomena now that wind velocities are increasingly more available. Estimating a relationship between SDEN and several climate drivers allows for the impact of climate variability on SDEN to be taken into account.