1. Introduction
Mountains are considered “water towers” since they provide water for both ecosystems and anthropogenic demands in downstream areas [1-3]. In Mediterranean mountains, this role becomes more important since melting water from the snowpack can constitute the main water resource for human consumption, agriculture, and hydroelectric production [
1]. The variability of the Mediterranean climate enhances the complexity of snow dynamics over these mountain regions [
2]. The mild winter temperatures in combination with the long dry sunny period undergo in: (i) a highly variable snowpack, in time and space, with the presence of several accumulation and melting cycles during the snow season [
3]; (ii) a shallow snowpacks with a characteristic patchy distribution [
4]; (iii) a high snowpack density [
5]; and (iv) a non-negligible evaposublimation flux from the snowpack to the atmosphere [
6]. This peculiar snow dynamics conditions the seasonality of the streamflow response in the head-water catchments [
2,
5]. Therefore, it is key to deeply understand accumulation and melting dynamics of the snowpack to be able to foresee changes in downstream streamflow and consequently, in the water resources availability [
7,
8]
Physics of snow dynamics have been widely studied [
9,
10,
11]. In a snowpack, during the accumulation phase snow water equivalent (SWE) increases, net energy inputs keep mainly negative and the average snowpack temperature decreases [
12,
13]. When these energy inputs become positive, due to external forcing, for instance radiation, snow turns into wet snow and the melting phase begins [
14,
15]. This melting period can be generally separated in three phases [
11]: (i) moistening; (ii) ripening ; and (iii) melting runoff. Traditionally, in-situ monitoring has been the principal tool used for observing and, therefore, understanding this snow dynamics [
16]. However, the intrinsic limitations of accessibility, measurement representativeness and area coverage makes nowadays remote sensing one of the techniques extensively used for snow monitoring [
10,
11,
17]. In addition, the possibility to use free satellite imagery on cloud platforms such as Google Earth Engine (GEE) is also expanding their use [
26,
27,
28]. Optical sensors, not just satellite sensor but also terrestrial photography, have been exploited to derive snow cover maps [
8,
13,
18,
19] and to compute snow albedo [
20,
21,
22], while, microwave information has been used to retrieve SWE from the snowpack [
23,
24,
25] , to estimate snow depth [
29], and to differentiate between dry and wet snow [
12,
25,
26].
In this last case, this differentiation is based on the fact that the presence of water in the snow cover generates high dielectric losses that provoke significant reduction in radar backscatter [
27]. This drop in the signal has been recently demonstrated that is considerably visible and sensible in the cross-polarization VH [
28,
29]. In particular, the new Sentinel-1 (S-1) SAR constellation, which has a dual polarization capability, very short revisit times, high spatial resolution, rapid product delivery and open access to all imagery [
30], has become widely used in different applications that use wet-snow information [
22,
23,
24,
25]. For instance, in [
31], the author used S-1 for characterizing wet snow to quantify mass balance over Greenland. In the same path, snowmelt over the ice sheets in Antarctica were characterized [
32]. In [
33], the authors used wet-snow identification for characterizing ice avalanches while in [
28], they have successfully characterized snowmelt phases using SAR in an Alpine environment. Moreover, in [
20], it is concluded that minimum SAR backscattering is reached at the end of the ripening phase and at the beginning of the melting-runoff phase. At that moment all liquid water content is held in the snowpack and starts to be released. From this onset, a monotonic increase in the backscattering is found until the snowpack is melted. Over these environments, snowpack is well consolidated during the winter, has thick depths and generally a unique snowmelt cycle during the springtime. Mediterranean mountains, located in semi-arid regions, are unlikely to have just one accumulation-melting cycle and suffer such slow transitions [
2,
4,
5,
34]. Indeed, in these areas the main characteristics are generally quicker snowpack changes, several melting cycles throughout the year and perhaps with the loss of some of the three defined melting phases due to this quick response [
3,
35].
From a hydrological point of view, the melting runoff onset provides the moment when melting water goes out of the snowpack and starts to contribute as liquid water into the system. Finding a relation between this onset and the streamflow response can help forecast streamflow behavior due to runoff melting [
36,
37]. In fact, the actual warming situation is altering this timing and increasing the need of a correct representation of this onset [
38,
39]. However, the paths followed by runoff melting water for reaching the stream are multiple and difficult to differentiate [
40]. In general, runoff melting water after leaving the snowpack arrives at the upper part of the soil [
41]. From there, meltwater infiltrates into the soil and/or accumulates to form a saturated zone at the base of the snowpack, moving toward a surface-water body following different paths that are conditioned by soil characteristics[
42]. If soil is unsaturated and water table is low, the water infiltrates and moves stream ward as subsurface flow. On the contrary, if soil is close to the saturation state, meltwater accumulates to form a basal saturated zone that develops within the snowpack through which water flows toward the stream [
11]. One can therefore conclude that melting dynamics are directly linked to baseflow response of the catchment.
In this context, the main objective of this study is twofold, first to explore the capability of C-band S-1 SAR imagery to capture multi seasonal wet snow dynamics over Mediterranean mountain areas and second, to apply a hydrological approach to assess the connection between this wet snow dynamics and streamflow response. In this work the main innovative parts are: i) the jointly use of terrestrial imagery and satellite information to deeply understand S-1 SAR backscattering response at the local scale; ii) the use of S-1 SAR imagery to characterize snow melting not only during the spring melting but throughout the hydrological year; and iii) the applicability of S-1 SAR imagery to identify runoff melting onset and to relate this behavior to changes in the streamflow dynamics at the catchment scale. The Sierra Nevada Mountain Range (southern Spain) is selected as pilot Mediterranean mountain area to carry out this study.
The paper is structured into five sections.
Section 2 presents the test sites at a local and catchment scales. It also includes the meteorological data available as well as the remote sensing information used. In section 3, we present the different steps of the proposed approach to detect wet snow conditions and the definition of melting cycles. The results are exposed in
Section 4, focusing on the local and catchment scale. In section 5, we discuss the main results. Finally, section 6 draws the conclusions of the work, by providing information on how the proposed methodology can be extended to other semi-arid mountains.
5. Discussion
The capability of S-1 SAR imagery to improve the understanding of wet snow dynamics over Mediterranean mountain areas was tested. We successfully connected melting dynamics with changes in the backscattering signal by applying the principles stated by [
28]. On top of that, we applied this approach not only to the final spring melting but also along all melting cycles that commonly appear throughout the year in these types of environments.
The characteristic of snow over these environments implied some technical challenges when adapting the general methodology to retrieve wet snow from SAR imagery. One of these challenges was the selection of a suitable reference image to remove common ground features that interfere in wet-snow detection. The lack of a long-lasting winter accumulation made difficult to find a full dry snow image, like commonly used in other studies [
61,
72]. Therefore, the selection of a dry summer image was chosen as reference. Several options were investigated following other authors [
83,
84]. Among them we observed that the soil dryness during summer month over our target catchments made almost no differences in the use of these different approaches.
Another issue was the definition of the threshold in the Δσ signal to discriminate wet snow. Here a value of -3.00 dB was chosen based on a wet-dry histogram analysis in the target area. This value agrees with other studies in alpine regions such as Austria or Argentina [
57,
72]. Another example of the use of a similar threshold is found in the Patagonian range in Argentina, a mountain closer to the Mediterranean than to the Alpine snow characteristics [
77]. These findings are also in line with the threshold defined by [
85,
86,
87], where wet snow was successfully analyzed in very different zones such as Finland or Canada. In addition, conditions such as the angle of incidence or the roughness of the snow also play an important role, in our case, angles of incidence of 35°-36° are in close proximity to the 30° analyzed by [
88], where the presence of snowpacks with high surface roughness were defined as an optimal range of values to define wet snow. In our case the shallow snowpack and the high snow roughness after some wind episodes visible in the terrestrial imagery support the selection of this threshold.
It is interesting to highlight that we have just used one of the two available orbits, which take place during the same day at different timing. Future analyses in the area should incorporate both orbits, being in this way able to capture hydrological differences between morning and evening snowpack. This daily changes are frequent in Mediterranean mountain like Sierra Nevada, especially during spring melting [
28].
The use of terrestrial photography also helped to deeply interpreted Δσ signal at the local scale. We were able to identify four types of LMs that triggered snowmelt-onset based on the classification proposed by [
3]. These LMs had different values depending on their timing throughout the year. This confirms the heterogeneity of the area, which can have different snowfalls events throughout the year, which ends up implying that we do not always have the minimum Δσ values in the long-lasting spring melting cycles.
This understanding of wet-snow dynamics helps to comprehend the melting-runoff onset maps at the catchment scale, which constitutes a useful tool for water managers to know the distribution of the melting dynamics. One limitation of these maps is that they are computed at the end of the snow season. They can be used to understand snow behavior but not to foresee changes in advance. On the contrary, the relationships found in section 4.2 goes a step further being able to be computed in quasireal real-time. On the one hand the depletion curves, with a common shape, defined a space where streamflow response and wet snow dynamics are linked using the maximum SWE of the analyzing melting cycle. Hence, knowing the maximum SWE achieved during the accumulation, a hypothetical curve can be drawn over this space to anticipate streamflow behavior. On the other hand, a mean delay of 21 days, between the beginning of the melting-runoff and streamflow peak was identified, narrowing down the water manager time to foresee runoff melting impact on the streamflow. More specifically, the proposed methodology can be applied in other sub-basins and be incorporated in operational hydrological modeling systems [
89,
90,
91].
Figure 1.
Location of Sierra Nevada within Spain; (b) Relief of the Sierra Nevada Mountain range with the location of the two study sites: Refugio Poqueira (yellow cross) and Poqueira Alto Catchment (blue polygon); (c) Zoom of the target areas.
Figure 1.
Location of Sierra Nevada within Spain; (b) Relief of the Sierra Nevada Mountain range with the location of the two study sites: Refugio Poqueira (yellow cross) and Poqueira Alto Catchment (blue polygon); (c) Zoom of the target areas.
Figure 2.
Workflow of the proposed approach where the main sections are identified in brackets. The inputs are SAR and MODIS images for section (i), the cycles are the results of the first step that is then used as an input to be interpreted at a local scale (ii) and at a sub-basin scale (iii).
Figure 2.
Workflow of the proposed approach where the main sections are identified in brackets. The inputs are SAR and MODIS images for section (i), the cycles are the results of the first step that is then used as an input to be interpreted at a local scale (ii) and at a sub-basin scale (iii).
Figure 3.
Scheme of melting cycle definition
Figure 3.
Scheme of melting cycle definition
Figure 5.
Definition of four types of local minima (red dots) identified for Mediterranean mountain area. VH represents the theoretical distribution of the dB of S-1 SAR backscattering of the cross-polarization VH. FSC represents the temporal distribution of snow cover fraction. Blue dots represent the maximum snow cover fraction. Periods without snow or with dry snow are also specified on the upper axis. The x axis represents the temporal variation over a hydrological year (Sep-Aug) of the different variables, the y axis represents the corresponding value of the dB of S-1 SAR backscattering (left y axis) and snow cover fraction (right y axis).
Figure 5.
Definition of four types of local minima (red dots) identified for Mediterranean mountain area. VH represents the theoretical distribution of the dB of S-1 SAR backscattering of the cross-polarization VH. FSC represents the temporal distribution of snow cover fraction. Blue dots represent the maximum snow cover fraction. Periods without snow or with dry snow are also specified on the upper axis. The x axis represents the temporal variation over a hydrological year (Sep-Aug) of the different variables, the y axis represents the corresponding value of the dB of S-1 SAR backscattering (left y axis) and snow cover fraction (right y axis).
Figure 6.
Temporal evolution of analyzed variables for the hydrological year 2018-2019 (a) Total area of snow pixels (grey line), dry snow pixels (blue line) and wet snow pixels (red line). Lines in the top and bottom of every panel identify an accumulation (blue) or melting (red) phase based on the MOD10A2 product. (b) temporal evolution of LM Summatory (blue line). (c) precipitation (blue bars) and snowfalls (gray bars) and the evolution of discharge (blue line), baseflow (dash blue line) and direct runoff (grey line) In all panels, the different melting cycles are highlighted by using a blue rectangle.
Figure 6.
Temporal evolution of analyzed variables for the hydrological year 2018-2019 (a) Total area of snow pixels (grey line), dry snow pixels (blue line) and wet snow pixels (red line). Lines in the top and bottom of every panel identify an accumulation (blue) or melting (red) phase based on the MOD10A2 product. (b) temporal evolution of LM Summatory (blue line). (c) precipitation (blue bars) and snowfalls (gray bars) and the evolution of discharge (blue line), baseflow (dash blue line) and direct runoff (grey line) In all panels, the different melting cycles are highlighted by using a blue rectangle.
Figure 7.
Temporal evolution of analyzed variables for the hydrological year 2018-2019 (a) Total area of snow pixels (grey line), dry snow pixels (blue line) and wet snow pixels (red line). Lines in the top and bottom of every panel identify an accumulation (blue) or melting (red) phase based on the MOD10A2 product. (b) the temporal evolution of the grouped temporally LM (blue line). (c) precipitation (blue bars) and snowfalls (gray bars) and the evolution of discharge (blue line), baseflow (dash blue line) and direct runoff (grey line) In all panels, the different melting cycles are highlighted by using a blue rectangle.
Figure 7.
Temporal evolution of analyzed variables for the hydrological year 2018-2019 (a) Total area of snow pixels (grey line), dry snow pixels (blue line) and wet snow pixels (red line). Lines in the top and bottom of every panel identify an accumulation (blue) or melting (red) phase based on the MOD10A2 product. (b) the temporal evolution of the grouped temporally LM (blue line). (c) precipitation (blue bars) and snowfalls (gray bars) and the evolution of discharge (blue line), baseflow (dash blue line) and direct runoff (grey line) In all panels, the different melting cycles are highlighted by using a blue rectangle.
Figure 8.
Temporal evolution of analyzed variables for the hydrological year 2018-2019 (a) Total area of snow pixels (grey line), dry snow pixels (blue line) and wet snow pixels (red line). Lines in the top and bottom of every panel identify an accumulation (blue) or melting (red) phase based on the MOD10A2 product. (b) the temporal evolution of the grouped temporally LM (blue line). (c) precipitation (blue bars) and snowfalls (gray bars) and the evolution of discharge (blue line), baseflow (dash blue line) and direct runoff (grey line) In all panels, the different melting cycles are highlighted by using a blue rectangle.
Figure 8.
Temporal evolution of analyzed variables for the hydrological year 2018-2019 (a) Total area of snow pixels (grey line), dry snow pixels (blue line) and wet snow pixels (red line). Lines in the top and bottom of every panel identify an accumulation (blue) or melting (red) phase based on the MOD10A2 product. (b) the temporal evolution of the grouped temporally LM (blue line). (c) precipitation (blue bars) and snowfalls (gray bars) and the evolution of discharge (blue line), baseflow (dash blue line) and direct runoff (grey line) In all panels, the different melting cycles are highlighted by using a blue rectangle.
Figure 9.
Depletion curves. the x-axis represents the dimensionless cumulative number of pixels contributing to melting runoff in each cycle (ΣLMo). The y-axis represents the dimensionless cumulative baseflow during each cycle. The 2016-2017 event is represented by the green lines. The 2017-2018 event is represented by the blue lines. The 2018-2019 events are represented by the purple line, the orange and brown line. The validation 2019-2020 event is represented by the red line. SWE max value for each cycle is also included in the legend. Red triangles represent the peak in baseflow during each event. Black dots represent the 25 % (●) of the total time of each cycle, black crosses represent 50% (x) of the total time of each cycle. and black plusses represent 75% (+) of the total time of each cycle.
Figure 9.
Depletion curves. the x-axis represents the dimensionless cumulative number of pixels contributing to melting runoff in each cycle (ΣLMo). The y-axis represents the dimensionless cumulative baseflow during each cycle. The 2016-2017 event is represented by the green lines. The 2017-2018 event is represented by the blue lines. The 2018-2019 events are represented by the purple line, the orange and brown line. The validation 2019-2020 event is represented by the red line. SWE max value for each cycle is also included in the legend. Red triangles represent the peak in baseflow during each event. Black dots represent the 25 % (●) of the total time of each cycle, black crosses represent 50% (x) of the total time of each cycle. and black plusses represent 75% (+) of the total time of each cycle.
Figure 10.
Baseflow and wet snow pixels temporal relationship. The x axis represents the number of wet snow pixels, the y axis represents the corresponding value of the base flow with its respective applied delay. (a) (b) (b).
Figure 10.
Baseflow and wet snow pixels temporal relationship. The x axis represents the number of wet snow pixels, the y axis represents the corresponding value of the base flow with its respective applied delay. (a) (b) (b).
Table 3.
Summary of the characteristics (beginning date, end date, cumulative precipitation, and cumulative snowfall) of snowmelt-accumulation events identified in each of the hydrological year analyzed.
Table 3.
Summary of the characteristics (beginning date, end date, cumulative precipitation, and cumulative snowfall) of snowmelt-accumulation events identified in each of the hydrological year analyzed.
Melting Cycle |
Beginning |
End |
Duration (days) |
Precipitation (mm) |
Snowfall (mm) |
Average daily Temperature (ºC) |
Average minimum daily Temperature (ºC) |
Mean SWE (mm) |
Maximum SWE (mm) |
2016- 2017 |
C-1 |
17-11-2016 |
05-12-2016 |
18 |
256 |
159 |
0.5 |
-3.0 |
85.2 |
140.2 |
C-2 |
11-12-2016 |
04-01-2017 |
24 |
92 |
67 |
1.3 |
-2.0 |
85.8 |
125.2 |
C-3 |
09-02-2017 |
04-05-2017 |
84 |
233 |
140 |
3.5 |
-0.6 |
23.9 |
78.7 |
2017- 2018 |
C-1 |
05-01-2018 |
23-01-2018 |
18 |
41 |
36 |
0.2 |
-4.0 |
14.3 |
26.4 |
C-2 |
29-01-2018 |
16-02-2018 |
18 |
32 |
31 |
2.9 |
-7.2 |
12.43 |
25.9 |
C-3 |
28-02-2018 |
09-08-2018 |
162 |
585 |
395 |
7.3 |
3.2 |
43.3 |
184.5 |
2018 – 2019 |
C-1 |
26-10-2018 |
07-11-2018 |
12 |
70 |
47 |
0.9 |
-3.5 |
14.4 |
34.5 |
C-2 |
13-11-2018 |
12-01-2019 |
60 |
139 |
85 |
3.7 |
0.4 |
13.8 |
52.3 |
C-3 |
30-01-2019 |
19-03-2019 |
48 |
85 |
72 |
2.5 |
-1.4 |
10.7 |
61.4 |
C-4 |
31-03-2019 |
23-06-2019 |
120 |
171 |
104 |
7 |
2.8 |
5.5 |
43.4 |