Snow cover as a morphogenic agent determining ground climate, landforms and runoff in the Valdecebollas massif, Cantabrian Mountains

Snowfalls are important meteorological events affecting the physical environment of the Cantabrian Mountains. This work analyzes the effects of snow on several elements such as relief, landforms, ground climate and snowmelt waters. The ground thermal regime and associated parameters were studied using temperature data loggers and satellite images and were described in combination with observed geomorphological processes and landforms. A geomorphological map was drawn up and trends in climate patterns and runoff were calculated. Ground temperature monitoring in warm years is not optimal, though allow to know the limit conditions for developing cold processes. Results show that geomorphological processes are not significant and that solifluction deriving from snowmelt, is the only active process in years without freeze or with thick snow cover. Snowfall evolution in recent decades in correlation with flow water and climate features provide the certainty that snow distribution also affects efficacy in runoff generation and moves the flow peak in rivers due to early snowmelt.

Another aspect of snow is its major role on mountain rivers water regime, e.g. in the Pyrenees (López-Moreno and García-Ruíz, 2004;López-Moreno et al. 2009;García-Ruíz et al. 2011;Sanmiguel et al., 2017).The dominant snow-fed fluvial regime is a characteristic of periglacial environments with a higher ratio between runoff and precipitation (Kane et al., 1998;Ballantyne, 2018).Snow accumulation produces a more efficient runoff regime that moves greater quantities of water and creates relevant geomorphological events (Woo, 2012).In the Duero/Douro catchment, which includes this massif, clear reductions in runoff and fluvial regime changes have already been recorded (Ceballos-Barbancho et al., 2008;Morán-Tejeda et al., 2010;Morán-Tejeda, 2012).Most of the research shows a very good correlation between decreases in snowfall and air temperature increases (Adam et al., 2009).Moreover, changes in the duration of snow cover also have consequences, bringing about alterations in microbial processes in soils (Edwards et al., 2007).
Nevertheless, most contributions to date on the Central Cantabrian Mountains have focused on the search for active processes (Castañón andFrochoso, 1994, 1998;Serrano and González-Trueba, 2004b;González-Trueba, 2007a;Santos-González et al., 2009;González-Trueba and Serrano, 2010b;Santos-González, 2011;Pellitero, 2013;Pisabarro et al., 2015Pisabarro et al., , 2017) ) and conclude that snow cover is the driver of freeze-thaw cycles (F/Tc) and seasonal frozen grounds (SFG).In this sense, the study of the ground thermal regime (GTR) facilitates the understanding of the capacity of the climate to trigger active periglacial processes and shows it to be a good indicator of snow cover characteristics (Zhang, 2005;French, 2007;Ballantyne, 2018), particularly in the Mediterranean mountains (Vieira et al., 2003;Andrés and Palacios, 2010;Oliva et al., 2014Oliva et al., , 2016) ) where the snow cover is highly variable in time because of the climatic characteristics of this region.
GTR studies in the Cantabrian Mountains confirmed the presence of an active periglacial belt above 2100 m (Oliva et al., 2016) according to the topoclimate.Nevertheless, this belt is beneath the threshold for mountain permafrost (French, 2007;Dobinski, 2011) except in Picos de Europa massif, where it may be found at the edge of ice patches (Serrano et al., 2011;Pisabarro et al., 2015) and in ice caves (Gómez Lende et al., 2014;Gómez-Lende, 2015).The other periglacial landforms were active in the Pleistocene and are good evidence of the past with permafrost e.g.rock glaciers in the Younger Dryas (YD) and Last Glacial Maximum (LGM) (Alonso, 1989;García de Celis, 1991;Serrano and Gutiérrez, 2000;Redondo-Vega et al., 2004, 2013;Santos-González et al., 2009;Pellitero et al., 2011;Rodríguez-Rodríguez et al., 2016;Oliva et al., 2018).This paper aims to synthesize the main processes of the physical environment in the Cantabrian Mountains and obtain new relationships among snow cover, ground climate, relief, landforms and hydric resources at a little studied site, the Pisuerga headwaters and Valdecebollas massif.

Study area
The Valdecebollas massif is in the North of Palencia at the source of the river Pisuerga, tributary of Duero catchment.The summit is at just 2143 m and is part of the western edge of the Mesozoic cover of the Cantabrian Mountains, touching Alto Campoo and Cantabria to the northeast (N 42º 58' 0,16"; W 4º 21' 55,94") (Fig. 1).This glacial landscape is appropriate for measuring these processes due to its topoclimatic features, glacial and periglacial landforms and recent cold processes such as ice waterfalls.The study area is north-facing over a glacial cirque with a tongue at the front that flowed to 1700 m in the LGM (Serrano et al., 2013).The glacial landforms are complex due to relief and geological structure.

P A P E R A C C E P T E D . P R E -P R I N T V E R S I O N
A wide variety of frontal and lateral moraines remain due to a succession of advances in three main stages, linked to a period since LGM to YD (Pellitero and Serrano, 2008;Serrano et al., 2013;Serrano and González-Trueba, 2004).The area shows clear evidence of fossil periglacial landforms above 1400 m as rock slopes in the Pisuerga valley (Pellitero and Serrano, 2008;Pisabarro, scheduled for 2019), a rock glacier in Alto Campoo (Serrano and Gutiérrez, 2000) Several asymmetric and dry NE-SW valleys are observed in Sierra Labra slopes (Nossin, 1959), splitting cryoplanation slopes in accordance with the definition by Ballantyne (2018).There are also covered ploughing boulders moved by a combination of solifluction and slumping (Nossin, 1959).Deep Pleistocene solifluction and gelifluction sheets were observed in 1.5 to 2 m profiles in stratified deposits in Tres Mares western face (Pisabarro, scheduled for 2019).
The rivers are fed by both snow and rainfall with flow peaks generally in spring due to the snowmelt, however the regime varies strongly due to the Mediterranean influence, which causes irregular warm episodes.During the study period (2015)(2016)(2017), mean winter air temperatures (December-March) of 2016-2017 were the warmest of the period 1961-2017 in this region increasing the historical average of winters (3,47ºC) in 2,02ºC in the closest long-term data (Fig. 2).Winter 2015-2016 was moderately warm (3,8ºC).Although this transitional climate has several warm and dry years, total precipitation is higher than 2,000 mm. and mean annual air temperature (MAAT) is 5,8ºC, minimum annual air temperature (Min AAT) is 2,8 ºC and maximum annual air temperature (Max AAT) is 9,1ºC (Fig. 3).
Logically, over several years, snowfalls are much more compensated and have a profound effect on the physical environment and socio-economic features of the area, e.g. the snowfalls of 2015 and 2018.Moreover, it is important to know that the area has been intensely affected by fires and overgrazing since the Middle Ages.

Methodology
The research is structured in three sections.Firstly, mean annual ground temperature (MAGT) features are analyzed, then active periglacial landforms are interpreted, and finally snow behaviour and runoff tendencies are compared to these of temperature and precipitation.

GTR features
GTR analysis was performed using thermal data loggers distributed in the Valdecebollas glacial cirque and peak (Table 1) buried at 0.1 m depth.The sensors used were iButton DS1921G with an accuracy of 0.5ºC and a range from -40 to +85ºC, which were programmed by software One Wire Viewer Maxim to take a reading with an acquisition frequency of 6 hours 6 hours over 11 consecutive months between December 2015 and November 2016, and between December 2016 and July 2017.The GTR provided the phases, which were obtained directly and show the evolution to provide a better understanding of thermic stabilities under snow cover or instabilities due to intense freeze or thaw (Delaloye, 2004;Zhang, 2005).Through these temperatures it was possible to know the number of F/Tc, which is the number of days with temperatures above and below 0ºC, the key to evaluating the SFG distribution (Frauenfeld et al., 2007).The freeze index (FI) (Fengqing and Yanwei, 2011) is obtained by adding daily temperatures below 0 ºC in each winter, 2015-2016 and 2016-2017.The depth of freeze was calculated using FI according to Washburn (1979); h(m)= +d; C L (latent heat 3.34•108 m 3 ); FI (freeze index °day); d (sensor depth 0,1 m) and K (conductivity values according to Eppelbaum et al. (2014)).The probability of a freeze ground is interpreted through the minimal temperatures period.If the temperature is below 0ºC, ± 0.5ºC, SFGs are certain and above > 0.5ºC unlikely.The MAGT map allows the identification of potential cold places for developing active periglacial landforms (Serrano et al., 2018).This was done using a simple interpolation of second grade in order to get a clear image of the MAGT distribution.The only reliable thermal image from the satellite Landsat was included with the thermal infrared obtained on 20th February 2017, in search of any signs of thermal anomalies.

Influence on the landforms
An inventory of periglacial landforms was made through fieldwork, photointerpretation and DEM high-resolution interpretation.These were then represented on a geomorphological map with a scale of 1 : 25,000.They were represented following the Swiss IGUL style (IGUL, 1996;Reynard and Lambiel, 2015) with some modifications, classified by colours depending on their origin.A similar method was used in the Fuentes Carrionas massif (Pellitero, 2014), whose activity was estimated according to GTR features, position and indicators of activity in the field such as absence of vegetation and lichens, signs of movement in the ground, geometry and alignments.The identification of the periglacial landforms was in accordance with French (2007) and Ballantyne (2018).Glacial landforms were also previously observed and compared (Hernández-Pacheco, 1944;Díaz-Martínez, 1993;Serrano and González-Trueba, 2004a;Pellitero and Serrano, 2008).
To see how snow affects runoff, a Spearman's non-parametric test was performed between the following two data series: (ii) The Pisuerga river headwater is the Valdecebollas glacial cirque and the closest gauging station with long series is at the Requejada reservoir (ETRS89: 375310; 4751674; 1083 m), located 15 km downstream.Q data correspond to the trends in the entrance flow.P and T data trends come from the meteorological stations located inside the catchment upstream of the reservoir (233 km 2 ) (Table 2).In addition, the historical Q, T, and P monthly distributions facilitate the construction of a linear regression with the series grouped by months and the extraction of Spearman's coefficient as a measure of positive or negative monthly trends in Q, P, and T over time.These trends permit us to interpret a change in the river flow regime (Beguería et al., 2003;López-Moreno et al., 2014).

GTR features
MAGT were extremely warm (Fig. 4), particularly in the 2015-2016 season.The lowest temperatures were recorded above 1980 m. Five sensors recorded days with temperatures below 0ºC.Sensor V5 had 15 days below 0ºC at 2134 m, and V6 had 22 days at 2,038 m in the winter of 2015-2016.In V1, V2, V3, V4 no ground temperatures (GT) below 0ºC were recorded because the altitude was less than 1910 m.Both sensors were oriented northwards (Fig. 4).The episode when most minimum temperatures occurred was in February 2016 during a long cold spell with -7ºC in the ground and a thin snow cover (Fig. 5).

P A P E R A C C E P T E D . P R E -P R I N T V E R S I O N
During 2016 temperatures did not fall below -2ºC.In the winter of 2016-2017 the threshold of -2ºC was crossed by V22 (7 days) and V66 (13 days) between December and January with GT of at most -8 ºC on the telenivometer.Between December 2016 and February 2017 GT were recorded below 0 ºC at sensor V55 in addition to the aforementioned V22 (14 days) and V66 (38 days).
With the exception of micro-cycles, which went undetected by the sensors, F/Tc were scarce, with 5 or less in whole data loggers with temperatures below 0ºC.All of these occurred together at the start of winter due to snow cover variability and successive changes in the atmosphere.The curtain effect was inferior to 5 months above 2000 m.FI is very low in 4 of the 5 mentioned data loggers (V5, V6, V22, V55), between 12 and 20.V66 reaches 80.7.Freeze depth, if was acting directly over the bedrock, would reach 50 cm in V66 due to the higher conductivity of conglomerates present on the top (Table 3).Data loggers V5, V6, V22 and V55 located above 1,982 m would reach 30 cm approximately.Below this threshold (1950 -2000 m) a clear incidence over the active layer can be ruled out.While it is true that temperatures in the winters analyzed were extremely warm, the thin snow cover favoured cold penetration.The possibility of SFG over very short periods is limited in the top parts (Table 3) depending on the snow cover, relief and cold intensity, e.g.January 2017 with ground temperatures below -2ºC (V55) and down to -6,5ºC (V66), and air temperatures around -8ºC in the telenivometer station at 1910 m.The GTR (Fig. 4) had only two clear phases, a warm regime adapted to air temperatures and a short equilibrium phase at around 0ºC in which the snow interferes with heat exchange (Ishikawa, 2003) called curtain effect (Outcalt et al., 1990).Transitional phases did not have a clear progressive trend.The curtain effect phase is short due to the limited snowfalls in these warm years.In V1, V2, V3, and V4 this phase was stable and close to zero reaching March but had positive MAGT in 2016.In V22 and V66 the curtain effect was observed from December 2016 to May 2017.It is feasible this change in GTR behavior indicates the moment of the snowmelt.V55 did not reach this month, the snow cover at the peak probably having disappeared earlier due to higher solar radiation or wind action.The satellite image with reflectance of thermal infrared band (Fig. 6) revealed the places with homogenous and thick snow cover.These reflected a colder response than irregular grounds, e.g.fossil landslide to the southwest of the image.Besides, the shadows of ridges clearly affect decreasing temperatures.The

P A P E R A C C E P T E D . P R E -P R I N T V E R S I O N
location of sensors with lowest GT, V66 and V55, reflects more energy, which means a thin snowpack because GT was above 0ºC on 20 th February 2017 (Fig. 4).

Influence on the landforms
The study area is full of fossil landforms from Pleistocene glacial activity, but also due to a small series of active landforms mostly linked to snow cover dynamics (Fig. 7).The large amounts of snow cover produce enough water melt to elevating moisture and mobilizing sediments along the slope.The active landforms with F/Tc and SFG action are limited to the area above 2000 m where the accumulation of snow is greater.Structural cryoturbations (Fig. 8) are clear above 2100 at the summit and under the north ridge.These processes trigger a certain organization of clasts, though not large enough to construct geometric structures (Ballantyne, 2018).Deposits are rich in clay with a high moisture level.Snow here is not stable due to the action of the wind.Below this summit area, where slopes increase between 2000 and 2100 m, terracettes have developed over shallow slope sediments, but only in the western slopes where the snowpack is reduced due to wind action.In the areas adjoining the terracettes there is brittle bedrock, mudrock lithology, erosion and the development of gullies with rills.There are also well developed solifluction lobes at one point favourable to cold at less than 2000 m.Other solifluction landforms are also visible in the top area.These can trigger the movement of ploughing boulders oriented downslope (Fig. 8), which are active around 2100-2050 m.These slope movements produced by snowmelt are created in many cases without F/Tc action (Matsuoka, 2001) , but they can also develop landforms deriving from snowmelt as sheets, bank-turf lobes, terraces with very ephemeral events of air freeze able to trigger surface forms such as needle ice above 1200 m.Sometimes freezing at night can develop pipkrake in combination with water melt moving slowly down the slopes.Without freeze features the water melt can also source enough moisture to trigger the movement of slopes and the reactiveness of certain asymmetrical valleys.Snow activity in the past has occurred with other nivation features that do not appear to be active.In the slopes above 1600 m it is possible to identify cryoplanation slopes, nivation hollows and snowbanks, abrasion surfaces and stratified deposits.Moreover, movements of disperse large blocks have taken place.In all of these landforms movement can be ruled out due to the arboreal and shrub vegetation.Landforms deriving from SFG are more disperse and only above slopes on which the snowpack does not accumulate.There are some block slopes and debris talus and cones.

Relationship between snow cover over runoff and river peak flow regimes
Days of snowfall at the pluviometry station of Stª Mª de Redondo (Fig. 9) fit well to the snow water equivalent modelled by Alonso-González et al. (2017) with a spatial resolution of 10 km and thickness recorded in situ by the telenivometer station.variable is able to approximate the main trends in volume of snow.In the period with overlapping there are similar trends so snowfall days data can be considered useful.
Snowfall days and the residuals of the linear regression of Q = aP+bT+c (Fig. 10) have a significant bivariate correlation at a correlation level of α < 0,01 (double tail) with a high Spearman's correlation coefficient ρ=0.412.The correlation is particularly strong after 1987 with fewer snowfalls.
On a monthly scale, the Spearman correlation coefficient indicates statistically significant changes in trends in water flow and temperature (Fig. 11).It exhibits significant decreasing trends in water flows between June and September.In the case of temperatures, trends differ greatly between winter and spring with respect to those of summer and early autumn.Nevertheless, there are no significant changes in monthly precipitations.

Discussion
MAGT recorded show a reduced incidence of cold that affects the ground to 30 or 50 cm in the coldest places during less than one month.Nevertheless, the depth of frozen grounds with the same air temperatures varies greatly depending on multiple factors rather than just the thickness and continuity of the snowpack, despite this being the most important variable (Washburn, 1979;Goodrich, 1982;Tanarro et al., 2010).SFG are likely above 2000 m and certain at 2100 m under these climatic features.This is in line with other mountains in the Iberian peninsula (Oliva et al., 2016) and 100 m higher than in the Picos de Europa between 2003 and 2007 and with Fuentes Carrionas between 2010 and 2012 (Pisabarro et al., 2017).
Snow cover is also the main agent, capable of modelling the slopes in the area, although its incidence is only relevant through the snowmelt.It saturates the ground moisture slope downwards and produces ephemeral solifluction on the slopes and nival creep after the snowmelt (Yamada, 1999).The huge amounts of snow in several years produce enough meltwater to bring about solifluction, even downslope from snow-patches (Kariya, 2002;Thorn and Hall, 2002) to 1200-1300 m.This is below the threshold without ephemeral frozen ground, at 1,980 m, although with sporadic needle ice or pipkrake phenomena (Fig. 8) to 1200-1300 m during seasonal atmospheric thermal inversion at the foot of the valleys.French (2007) limits solifluction to 50 cm of the active layer only under F/Tc and slope movements between 1 and 10 cm/yr.This depth is similar to the activity around active periglacial landforms found depending on F/Tc, such as terracettes and cryoturbations.These landforms only remain active over 2000 m, as do terracettes in the Central Range (de Marcos et al., 2017).Solifluction lobes located at 1932 m involve thicker active ground, as a result of additional snow melt water inputs (Fig. 8).
Since snow cover protects the ground and landforms linked to F/Tc are scarce, e.g.frost weathering, frost shattering, debris fall.The development of gullies, rills and naked ground (Figs. 6 and  7) in the ridges and nearby slopes with unconsolidated lithologies is caused by the erosive intensity of cold and the push movement of snow cover capable of breaking the rocks (Carrera and Valcárcel, 2010).Even the periglacial deposits increase the nival action (Palacios and García-Sánchez, 1997) and several nival hollows may remain active (García-Sancho et al., 2001).Thermal reflectance (Fig. 6) shows evidence of a lower temperature on the surface if the snowpack is homogeneous, unlike the rough ground in the southwest of the image.Here there is a likely landslide without clear chronology and modified by solifluction movements with features of active-layer failure under Pleistocene permafrost deglaciation conditions.These conditions seem likely if we consider the actual warm conditions in 2017 MAGT of 0.81ºC in sensor V66, close to 0ºC (French, 2007;Ballantyne, 2018).
Although 'snowfall days' is an inexact variable, the correlation is strong enough to suggest that snowfalls are determinant in producing runoff, as occurs in Pyrenees (López-Moreno and García-Ruiz, P A P E R A C C E P T E D .P R E -P R I N T V E R S I O N 2004).While it is true that there are doubts as to whether pluviometers underestimate precipitations in snowfall (Buisán et al., 2017), from the data the influence of the snowpack over runoff is clear.Nor is the karstic environment at some points responsible for water loss by catchment changes (Alcalde, 1981) in the Pisuerga headwaters where the Cobre cave drives water to the Pisuerga.On a monthly scale, P, Q and T, significant changes can be seen as Morán-Tejeda (2012) stated in this same reservoir using management data.The Spearman's coefficient of Q decreased greatly in spring and early summer and increased in winter.Q change in spring may be slightly related with P change, but the influence of T seems to be determinant since its decrease is more marked and moves the snowmelt flow peak forward due to the snow melt.This behaviour is in common with other catchments in the Duero/Douro system (Ceballos- Barbancho et al., 2008;Morán-Tejeda et al., 2010).This may be evidence of a less important role of snow in the river regime, with an earlier flow peak.Evidence of this has been found in Pyrennean rivers (García-Ruíz et al., 2015) and could develop less effective geomorphologic activity in spring, as occurs in the Central Range (Palacios et al., 2003) if we take into account the irregularity of snowfalls, which in some years produced exceptional thicknesses (2014-15 and 2017-18).

Conclusions
Snow cover is the main variable determining GTR, active periglacial landforms and hydric resources in the Valdecebollas massif and similar behaviors could be expected on areas with similar characteristics (i.e.Cantabrian Mountains).
Even in warm years with a scarce number of snowfalls and snowpack permanence this can limit frost action, and the only active periglacial landforms linked to GTR are above 2000 m in places unfavourable to snowpack accumulation and continuity or in dry and cold years.
Hydric resources are also determined by the snowpack, which dictates the efficacy of solid precipitations by generating more or less runoff.By means of the monthly evolution of flows we can also estimate that the snowpack melts earlier and earlier and can modify peak flows.

Figure 3 .
Figure 3. Mean climate features at Valdecebollas telenivometer (precipitation estimated by snow wáter equivalent values.Rainfalls are excluded, so real values are greater. (i) Number of days with snow precipitation between 1961 and 2005.P A P E R A C C E P T E D .P R E -P R I N T V E R S I O N (ii) Residual values of linear regression among water flow (Q), precipitation (P) and air temperatures (T) in the Pisuerga catchment.These residuals are the values of Q not explained by P and T.Data used to construct the above series:(i) The pluviometry station of St. M. de Redondo (ETRS89: 382986; 4760666; 1200 m) located 6 km away continuously acquiring observations between 1961 and 2005.It is the only snowfall data series close to the area before 2008.The snowfall days were observed and collected by a person As this is not a very reliable variable, it was checked for trends using the Factorial SnowModel (1980Model ( -2014) )  (Alonso-González et al., 2017) and with the snow water equivalent values recorded by the Valdecebollas telenivometer (ETRS89: 389178; 4758644; 1910 m), which has been active since 2008 thanks to the Ebro Catchment Authorities, although this area belongs to the Duero Catchment Authorities.

Figure 5 .
Figure 5. MAGT in the study area.

Figure 6 .
Figure 6.Left: multispectral image.All the image is covered by snow.Right: classification of the reflectance of the spectral band of thermal infrared (B-10).Satellite image.Landsat 8. 23 rd February 2017.

Figure 7 .
Figure 7. Geomorphological map of Valdecebollas massif.P A P E R A C C E P T E D .P R E -P R I N T V E R S I O N Figure 9. Standardized comparison of different snow variables with the objective of checking if snowfall daysvariable is able to approximate the main trends in volume of snow.In the period with overlapping there are similar trends so snowfall days data can be considered useful.

Figure 10 .
Figure 10.Temporal evolution of residuals of linear regression (Q = aP+bT+c) and snowfall days in Stª Mª de Redondo.P A P E R A C C E P T E D .P R E -P R I N T V E R S I O N

Table 1 .
Characteristics of the sites with temperature data loggers installeg.

Table 3 .
Ground thermal results synthesis.