Snow sensitivity to temperature and precipitation change during compound cold–hot and wet–dry seasons in the Pyrenees

The Mediterranean Basin has experienced one of the highest warming rates on earth during the last few decades, and climate projections predict water scarcity in the future. Mid-latitude Mediterranean mountain areas, such as the Pyrenees, play a key role in the hydrological resources for the highly populated lowland areas. However, there are still large uncertainties about the impact of climate change on snowpack in the high mountain ranges of this region. Here, we perform a snow sensitivity to temperature and precipitation change analysis of the Pyrenean snowpack (1980–2019 period) using five key snow–climatological indicators. We analyzed snow sensitivity to temperature and precipitation during four different compound weather conditions (cold–dry (CD), cold–wet (CW), warm–dry (WD), and warm–wet (WW)) at low elevations (1500 m), mid elevations (1800 m), and high elevations (2400 m) in the Pyrenees. In particular, we forced a physically based energy and mass balance snow model (FSM2), with validation by ground-truth data, and applied this model to the entire range, with forcing of perturbed reanalysis climate data for the period 1980 to 2019 as the baseline. The FSM2 model results successfully reproduced the observed snow depth (HS) values (R²>0.8), with relative root mean square error and mean absolute error values less than 10 % of the observed HS values. Overall, the snow sensitivity to temperature and precipitation change decreased with elevation and increased towards the eastern Pyrenees. When the temperature increased progressively at 1 ^∘C intervals, the largest seasonal HS decreases from the baseline were at +1 ^∘C. A 10 % increase in precipitation counterbalanced the temperature increases (≤1 ^∘C) at high elevations during the coldest months because temperature was far from the isothermal 0 ^∘C conditions. The maximal seasonal HS and peak HS max reductions were during WW seasons, and the minimal reductions were during CD seasons. During WW (CD) seasons, the seasonal HS decline per degree Celsius was 37 % (28 %) at low elevations, 34 % (30 %) at mid elevations, and 27 % (22 %) at high elevations. Further, the peak HS date was on average anticipated for 2, 3, and 8 d at low, mid, and high elevation, respectively. Results suggest snow sensitivity to temperature and precipitation change will be similar at other mid-latitude mountain areas, where snowpack reductions will have major consequences for the nearby ecological and socioeconomic systems.

1 Introduction

Snow is a key element of the earth's climate system (Armstrong and Brun, 1998) because it cools the planet (Serreze and Barry, 2011) by altering the surface energy balance (SEB), increasing the albedo, and modulating surface and air temperatures (Hall, 2004). Northern hemispheric snowpack patterns have changed rapidly in recent decades (Hammond et al., 2018; Hock et al., 2019; Notarnicola et al., 2020). It is crucial to improve our understanding of the timing of snow ablation and snow accumulation due to changing climate conditions because snowpack affects many nearby social and environmental systems. From a hydrological point of view, snowmelt controls the mountain runoff rate during the spring (Barnett et al., 2005; Adam et al., 2009; Stahl et al., 2010), river flow magnitude and timing (Morán-Tejeda et al., 2014; Sanmiguel-Vallelado et al., 2017), water infiltration and groundwater storage (Gribovszki et al., 2010; Evans et al., 2018), and transpiration rate (Cooper et al., 2020). The presence and duration of snowpack affect terrestrial ecosystem dynamics because the snow ablation date affects photosynthesis (Woelber et al., 2018), forest productivity (Barnard et al., 2018), freezing and thawing of the soil (Luetschg et al., 2008; Oliva et al., 2014), and thickness of the active layer in permafrost environments (Hrbáček et al., 2016; Magnin et al., 2017). Snowpack also has remarkable economic impacts. For example, the snowpack at high elevations and in surrounding areas determines the economic success of many mountain ski resorts (Scott et al., 2003; Pons et al., 2015; Gilaberte-Búrdalo et al., 2017). Changes in the snowpack of mountainous regions also influence associated lowland areas because it affects the availability of snow meltwater that is used for water reservoirs, hydropower generation, agriculture, industries, and other applications (e.g., Sturm et al., 2017; Beniston et al., 2018).

Mid-latitude snowpacks have among the highest snow sensitivities worldwide (Brown and Mote, 2009; López-Moreno et al., 2017, 2020b). In regions at high latitudes or high elevations, increasing precipitation can partly counterbalance the effect of increases in temperature on snow cover duration (Brown and Mote, 2009). Climate warming decreases the maximum and seasonal snow depth (HS), snow water equivalent (SWE) (Trujillo and Molotch, 2014; Alonso-González et al., 2020b; López-Moreno et al., 2013a, 2017), and fraction total precipitation as snowfall (snowfall ratio; e.g., Mote et al., 2005; Lynn et al., 2020; Jennings and Molotch, 2020; Marshall et al., 2019), and it also delays the snow onset date (Beniston, 2009; Klein et al., 2016). However, warming can slow the early snow ablation rate in the season (Pomeroy et al., 2015; Rasouli et al., 2015; Jennings and Molotch, 2020; Bonsoms et al., 2022; Sanmiguel-Vallelado et al., 2022) because of the earlier HS and SWE peak dates (Alonso-González et al., 2022), which coincide with periods of low solar radiation (Pomeroy et al., 2015; Musselman et al., 2017a).

The Mediterranean Basin is a region that is critically affected by climate change (Giorgi, 2006), being densely populated (>500 million inhabitants) and affected by intense anthropogenic activity. Warming of the Mediterranean Basin will accelerate for the next few decades, and temperatures will continue to increase in this region during the warm months (Knutti and Sedlacek, 2013; Lionello and Scarascia, 2018; Cramer et al., 2018; Evin et al., 2021; Cos et al., 2022), increasing atmospheric evaporative demands (Vicente-Serrano et al., 2020) and drought severity (Tramblay et al., 2020), leading to water scarcity over most of this region (García-Ruiz et al., 2011). Mediterranean mid-latitude mountains, such as the Pyrenees, where this research is focused, are the main runoff-generation zones of the downstream areas (Viviroli and Weingartner, 2004) and provide most of the water used by major cities in the lowlands (Morán-Tejeda et al., 2014).

Snow patterns in the Pyrenees have high spatial diversity (Alonso-González et al., 2019) due to internal climate variability in mid-latitude precipitation (Hawkins and Sutton, 2010; Deser et al., 2012) and high interannual and decadal variability in precipitation in the Iberian Peninsula (Esteban-Parra et al., 1998; Peña-Angulo et al., 2020), as well as the abrupt topography and the different mountain exposure to the Atlantic air masses (Bonsoms et al., 2021a). Thus, snow accumulation per season is almost twice as much in the northern slopes as in the southern slopes (Navarro-Serrano and López-Moreno, 2017), and there is a high interannual variability of snow in regions at lower elevations (Alonso-González et al., 2020b) and in the southern and eastern regions of the Pyrenees (Salvador Franch et al., 2014; Salvador-Franch et al., 2016; Bonsoms et al., 2021b). Since the 1980s, the energy available for snow ablation has significantly increased in the Pyrenees (Bonsoms et al., 2022), and winter snow days and snow accumulation have non-statically significantly increased (Buisan et al., 2016; Serrano-Notivoli et al., 2018; López-Moreno et al., 2020a; Bonsoms et al., 2021a) due to the increasing frequency of positive west and southwest advections (Buisan et al., 2016). The 21st century climate projections for the Pyrenees anticipate a temperature increase of more than 1 to 4 ^∘C (relative to 1986–2005) and an increase (decrease) in precipitation by about 10 % for the eastern (western) regions during winter and spring (Amblar-Francés et al., 2020). Therefore, changes in snow patterns in regions with high elevations are uncertain because winter snow accumulation is affected by precipitation (López-Moreno et al., 2008), and Mediterranean Basin winter precipitation projections have uncertainties of up to 80 % of the total variance (Evin et al., 2021).

Previous studies in the central Pyrenees (López-Moreno et al., 2013a), Iberian Peninsula mountain ranges (Alonso-González et al., 2020b), and mountain areas that have Mediterranean climates (López-Moreno et al., 2017) demonstrated that snowpack sensitivity to changes in climate is mostly controlled by elevation. Despite the impact of climate warming on mountain hydrological processes, there is limited understanding of snow sensitivity to temperature and precipitation changes and seasonality of mid-latitude Mediterranean mountain snowpacks. Some studies have reported different snowpack sensitivities during wet and dry years (López-Moreno et al., 2017; Musselman et al., 2017b; Rasouli et al., 2022; Roche et al., 2018). However, the sensitivity of snow during periods when there are seasonal compound weather (temperature and precipitation) conditions has not been analyzed yet. The high interannual variability of the Pyrenean snowpack, which is expected to increase according to climate projections (López-Moreno et al., 2008), indicates a need to examine snowpack sensitivity to temperature and precipitation change focusing on the year-to-year variability. Warm seasons in the Mediterranean Basin require special attention because these are likely to increase in the future (e.g., Vogel et al., 2021; De Luca et al., 2020; Meng et al., 2022). Further, the occurrence of different HS trends at mid- and high-elevation areas of this range (López-Moreno et al., 2020a) suggests that elevation and spatial factors contribute to the wide variations in the sensitivity of snow to the climate.

Therefore, the main objective of this research is to quantify snow (accumulation, ablation, and timing) sensitivity to temperature and precipitation change during compound temperature and precipitation seasons in the Pyrenees.

2 Geographical area and climate setting

The Pyrenees is a mountain range located in the north of the Iberian Peninsula (southern Europe; 42–43^∘ N to 2^∘ W–3^∘ E) that is aligned east to west between the Atlantic Ocean and the Mediterranean Sea. The highest elevations are in the central region (Aneto, 3404 m), and elevations decrease towards the west and east (Fig. 1). The Mediterranean Basin, including the Pyrenees, is in a transition area and is influenced by the continental climate and the subtropical temperate climate. Precipitation is mostly driven by large-scale circulation patterns (Zappa et al., 2015), the jet-stream oscillation during winter (Hurrell, 1995), and land–sea temperature differences (Tuel and Eltahir, 2020). During the summer, the northward movement of the Azores high-pressure region brings stable weather, and precipitation is mainly convective at that time (Xercavins, 1985). Precipitation is highly variable depending on mountain exposure to the main circulation weather types; it ranges from about 1000 mm yr⁻¹ to about 2000 mm yr⁻¹ (in the mountain summits), with lower levels in the east and south (Cuadrat et al., 2007). There is a slight disconnection of the general climate circulation towards the eastern Pyrenees, where the Mediterranean climate and East Atlantic–West Russia (EA–WR) oscillations have greater effects on snow accumulation (Bonsoms et al., 2021a). In the southern, western, and central massifs of the range, the Atlantic climate and the negative North Atlantic Oscillation (NAO) phases regulate snow accumulation (W and SW wet airflows; López-Moreno, 2005; López-Moreno and Vicente-Serrano, 2007; Buisan et al., 2016; Alonso-González et al., 2020a). In the northern slopes, the positive phases of the Western Mediterranean Oscillation (WeMO) linked with NW and N advections trigger the most episodes of snow accumulation (Bonsoms et al., 2021a). The seasonal snow accumulation in the northern slopes is almost double the amount (about 500 cm more) as in the southern slopes at an elevation of about 2000 m (Bonsoms et al., 2021a). The temperature/elevations gradient is about 0.55 ^∘C per 100 m (Navarro-Serrano and López-Moreno, 2017), and the annual 0 ^∘C isotherm is at about 2750 to 2950 m (López-Moreno and García-Ruiz, 2004). Net radiation and latent heat flux govern the energy available for snow ablation; the former heat flux increases at high elevations and the latter towards the east (Bonsoms et al., 2022).

Figure 1 (a) Study area. Pyrenean massifs grouped by sectors for (b) low, (c) mid, and (d) high elevation. The white dots indicate the locations of the automatic weather stations (AWSs) shown in Table 1. Massif delimitation is based on the spatial regionalization of the Système d’Analyse Fournissant des Renseignements Adaptés à la Nivologie (SAFRAN) system, which groups massifs according to topographical and meteorological characteristics (modified from Durand et al., 1999).

3 Data and methods

3.1 Snow model

Snowpack was modeled using a physical-based snow model, the Flexible Snow Model (FSM2; Essery, 2015). This model resolves the SEB and mass balance to simulate the state of the snowpack. FSM2 is an open-access model and is available at https://github.com/RichardEssery/FSM2 (last access: 16 December 2022). Previous studies have tested FSM2 (Krinner et al., 2018) and its application in different forest environments (Mazzotti et al., 2021) and hydroclimatological mountain zones such the Andes (Urrutia et al., 2019), the Alps (Mazzotti et al., 2020), Colorado (Smyth et al., 2022), the Himalayas (Pritchard et al., 2020), Iberian Peninsula mountains (Alonso-González et al., 2020b, 2022), and Lebanese mountains (Alonso-González et al., 2021), providing confident results. FSM2 requires forcing data of precipitation, air temperature, relative humidity, surface atmospheric pressure, wind speed, incoming shortwave radiation (SW_inc), and incoming longwave radiation (LW_inc). We have evaluated different FSM2 model configurations (not shown) without remarkable differences in the accuracy and performance metrics. Thus, the FSM2 configuration included in this work estimates snow cover fraction based on a linear function of HS and albedo based on a prognostic function, with increases due to snowfall and decreases due to snow age. Atmospheric stability is calculated as a function of the Richardson number. Snow density is calculated as a function of viscous compaction by overburden and thermal metamorphism. Snow hydrology is estimated by gravitational drainage, including internal snowpack processes, runoff, refreeze rates, and thermal conductivity. Table S1 summarizes the FSM2 configuration and the FSM2 compile numbers.

3.2 Snow model validation

FSM2 configuration was validated by in situ snow records of four automatic weather stations (AWSs) that were at high elevations in the Pyrenees. Precipitation in mountainous and windy regions is usually affected by undercatch (Kochendorfer et al., 2020). Thus, the instrumental records of precipitation were corrected for undercatch by applying an empirical equation validated for the Pyrenees (Buisán Sanz et al., 2019). Precipitation type was classified by a threshold method (Musselman et al., 2017b; Corripio and López-Moreno, 2017): snow when the air temperature was below 1 ^∘C and rain when the air temperature was above 1 ^∘C, according to previous research in the study area (Corripio and López-Moreno, 2017). The LW_inc heat flux of the AWSs (Table 1) was estimated as previously described (Corripio and López-Moreno, 2017). Due to the wide instrumental data coverage (99.3 % of the total dataset), gap filling was not performed. The HS records were measured at each 30 min interval using an ultrasonic snow depth sensor. The meteorological data used in the validation process were provided and are managed by the local Meteorological Service of Catalonia (https://www.meteo.cat/wpweb/serveis/formularis/peticio-dinformes-i-dades-meteorologiques/peticio-de-dades-meteorologiques/; last access: 14 January 2021). Quality checking of the data was performed using an automatic error filtering process in combination with a climatological, spatial, and internal coherency control defined by the Servei Meteoròlogic de Catalunya (2011).

Table 1 Characteristics of the four AWSs.

Model accuracy was estimated based on the mean absolute error (MAE) and the root mean square error (RMSE), and model performance was estimated by the coefficient of determination (R²). The MAE and the RMSE indicate the mean differences of the modeled and observed values.

3.3 Atmospheric forcing data

We forced FSM2 with the open-access climate reanalysis dataset provided by Vernay et al. (2022), which consists of the modeled values from the SAFRAN meteorological analysis. FSM2 was run at an hourly resolution for each massif, each elevation range, and each climate baseline and perturbation scenario from 1980 to 2019. The SAFRAN system provides data for homogeneous meteorological and topographical mountain massifs every 300 m from 0 to 3600 m (Durand et al., 1999; Vernay et al., 2022). We analyzed three elevation bands: low (1500 m), middle (1800 m), and high (2400 m). Precipitation type was classified using the same threshold approach used for model validation. Atmospheric emissivity was derived from the SAFRAN LW_inc and air temperature. SAFRAN was forced using numerical weather prediction models (ERA-40 reanalysis data from 1958 to 2002 and Action de Recherche Petite Échelle, Grande Échelle (ARPEGE) from 2002 to 2020). Meteorological data were calibrated, homogenized, and improved by in situ meteorological observation data assimilation (Vernay et al., 2022). Durand et al. (1999, 2009a, b) provided further technical details of the SAFRAN system. Previous studies have used the SAFRAN system for long-term HS trends (López-Moreno et al., 2020a), extreme snowfall (Le Roux et al., 2021), and snow ablation analysis (Bonsoms et al., 2022). The SAFRAN system has been extensively validated for the meteorological modeling of continental Spain (Quintana-Seguí et al., 2017), France (Vidal et al., 2010), and alpine snowpack climate projections (Verfaille et al., 2018) among other works.

3.4 Snow sensitivity to temperature and precipitation change analysis

Snow sensitivity to temperature and precipitation change was analyzed using a delta change methodology (López-Moreno et al., 2008; Beniston et al., 2016; Musselman et al., 2017b; Marty et al., 2017; Alonso-González et al., 2020b; Sanmiguel-Vallelado et al., 2022). In this method, air temperature and precipitation were perturbed for each massif and elevation range based on the historical period (1980–2019). Air temperature was increased from 1 to 4 ^∘C at 1 ^∘C intervals, assuming an increase of LW_inc accordingly (Jennings and Molotch, 2020). Precipitation was changed from −10 % to +10 % at 10 % intervals, in accordance with climate model uncertainties and the maximum and minimum precipitation projections for the Pyrenees (Amblar-Francés et al., 2020).

3.5 Snow climate indicators

Snow sensitivity to temperature and precipitation change was analyzed using five key indicators: (i) seasonal average HS, (ii) seasonal maximum absolute HS peak (peak HS max), (iii) the date of the maximum HS (peak HS date), (iv) the number of days with HS >1 cm on the ground (snow duration), and (v) daily average snow ablation per season (snow ablation, hereafter). Snow ablation was calculated as the difference between the maximum daily HS recorded on 2 consecutive days (Musselman et al., 2017a), and only days with decreases of 1 cm or more were recorded. Some seasons had more than one peak HS; for this reason, the peak HS date was determined after applying a moving average of 5 d. All indicators were computed according to massif and elevation range.

3.6 Definitions of compound temperature and precipitation seasons

The snow season was from 1 October to 30 June (inclusive). Snow duration was defined by snow onset and snow ablation date in situ observations (Bonsoms, 2021a) and results from the baseline scenario snow duration presented in this work. A “compound temperature and precipitation season” (season type) was assessed based on each massif and the elevation historical climate record (1980–2019) using a joint-quantile approach (Beniston and Goyette, 2007; Beniston, 2009; López-Moreno et al., 2011b). Compound season types were defined according to López-Moreno et al. (2011b) based on the seasonal 40th percentiles (T40 for temperature and P40 for precipitation) and the seasonal 60th percentiles (T60 and P60). There were four types of seasons based on seasonal temperature (Tseason) and seasonal precipitation (Pseason) data:

• cold and dry (CD), Tseason ≤ T40 and Pseason ≤ P40;

• cold and wet (CW), Tseason ≤ T40 and Pseason ≥ P60;

• warm and dry (WD), Tseason > T40 and Pseason ≤ P40;

• warm and wet (WW), Tseason > T60 and Pseason > P60.

All remaining seasons were classified as having average (Avg) temperature and precipitation. Note that the number of compound season types is different depending on the Pyrenees massif (Fig. S1). However, by applying the joint-quantile approach described above, we are comparing the snow sensitivity to temperature and precipitation change between similar climate conditions, independently of where each compound season type was recorded.

3.7 Spatial regionalization

We have examined spatial differences in the snow sensitivity to temperature and precipitation change by compound season types. Following previous studies, massifs were grouped into four sectors by applying a principal component analysis (PCA) (i.e., López-Moreno et al., 2020b; Matiu et al., 2020, among others). We applied a PCA over HS data for each month, year, massif, and elevation. Massifs were grouped into fours sectors depending on the maximum correlation with PC1 and PC2 scores (see Fig. S2). The number of season types per sector are shown in Fig. S3, and the spatial regionalization is presented in Fig. 1.

4 Results

We validate FSM2 in Sect. 4.1. Subsequently, we analyze the snow sensitivity to temperature and precipitation change based on five snow climate indicators, namely the seasonal HS, peak HS max, peak HS date, snow duration, and snow ablation. Compound season types show similar relative importance of the snow sensitivity to temperature and precipitation change regardless of the Pyrenean sector. For this reason, our results have been focused on seasonal snow changes due to increments of temperature, elevation, and compound season type. These are the key factors that ruled the snow sensitivity to temperature and precipitation change, and an accurate analysis is provided in Sect. 4.2. Spatial differences of the snow sensitivity to temperature and precipitation change during compound season types are examined in Sect. 4.3.

4.1 Snow model validation

Our snow model validation analysis (Figs. 2 and 3) confirmed that FSM2 accurately reproduces the observed HS values. On average, FSM2 had a R² greater than 0.83 for all stations. In general, the snow model slightly overestimated the maximum HS values. The highest R² values were at A4 and A2 (R²=0.85 in both stations), and the lowest were at A3 and A1 (R²=0.79 and R²=0.82, respectively). The highest accuracy was at A4 (RMSE =18.5 cm, MAE =8.9 cm), and the largest errors were at A2 (RMSE =45.8 cm, MAE =29.0 cm).

Figure 2 Time series of the observed (gray) and modeled (black) HS values at the four AWSs.

Figure 3 Regression analysis of observed (x axis) and simulated (y axis) HS values.

4.2 Snow sensitivity to temperature and precipitation change

We then determined seasonal HS profiles for each perturbed climate scenario and compound season type (Fig. 4). The results show a non-linear response between seasonal HS loss and temperature increase. When the temperature increased at 1 ^∘C intervals, the largest relative seasonal HS decrease from the baseline was at +1 ^∘C for all elevations and all compound season types. High-elevation areas had lower seasonal HS variability between compound season types than low elevations (Fig. S4). At low elevations, snow was greater during CW seasons than other seasons. All the snowpack-perturbed scenarios indicated that snowpack decreased for all elevations under warming climate scenarios. Snowpack sensitivity to temperature and precipitation change depended on the compound season type (Figs. 5 and 6). At low elevations, the seasonal changes in HS ranged from −37 % ^∘C⁻¹ (WW) to −28 % ^∘C⁻¹ (CD). For mid-elevation ranges, there were no remarkable differences among compound season types (Table 2), and the seasonal HS changes ranged from −34 % (WW) to −30 % (CW) per degree Celsius increase. Low and mid elevations had greater snowpack reductions than high elevations. In the latter, a 10 % increase in precipitation counterbalanced a temperature increase of about 1 ^∘C, and there were no remarkable differences in the seasonal HS from the baseline scenario, especially in the coldest months of the season (Figs. S5 and S6). The maximum seasonal HS sensitivity to temperature and precipitation change was during WD seasons (27 % ^∘C⁻¹) and the minimum was during CW seasons (−22 % ^∘C⁻¹).

Figure 4 Anomalies of seasonal HS for low, mid, and high elevation (rows); compound season type (columns); and different temperature increases (colors).

Table 2 Average and seasonal HS and peak HS sensitivity to temperature and precipitation change during the four different compound temperature and precipitation seasons at three different elevations.

At low and mid elevations, the peak HS max was greatest during WW seasons (−24 % ^∘C⁻¹) and lowest during the CD and WD seasons (−17 % ^∘C⁻¹ for both). At high elevations, there were no clear differences in the peak HS max for the different seasons. The maximum peak HS max was during WD seasons (−16 % ^∘C⁻¹) and the minimum was during CD seasons (−14 % ^∘C⁻¹).

We also determined the average seasonal snow duration for each elevation range and the compound season type for different temperature increases (Table 3 and Fig. 5c). The minimum snow duration was during CW seasons (−13 % ^∘C⁻¹ at low elevations, −10 % ^∘C⁻¹ at mid elevations, −5 % ^∘C⁻¹ at high elevations). At low elevations, the snow duration was the most sensitive during WW seasons (−17 % ^∘C⁻¹). On the contrary, at mid and high elevations, the snow duration was the most sensitive during WD seasons (−13 % ^∘C⁻¹ at mid elevations and −8 % ^∘C⁻¹ at high elevations).

Figure 5 Anomalies of seasonal HS (a), peak HS max (b), and snow duration (c) for different temperature increases relative to the baseline at three different elevations during the four different compound season types. The solid black lines within each boxplot are the average. Lower and upper hinges correspond to the 25th and 75th percentiles, respectively. The whisker is a horizontal line at a 1.5 interquartile range of the upper and lower quartile. Dots represent the outliers. Data are grouped by season, compound season type, increment of temperature, precipitation variation, elevation, and massif.

The peak HS date occurred earlier due to warming, independently of precipitation changes. During WD seasons, the peak HS date per degree Celsius was anticipated by 3 d per season at low elevations, 3 d at mid elevations, and 6 d at high elevations; during CD seasons, the peak HS date per degree Celsius was anticipated by 4 d per season at low elevations, 5 d at mid elevations, and 9 d at high elevations. In low- and mid-elevation areas, if the temperature increase was no more than about 1 ^∘C above the baseline, there was little change in the peak HS date (Fig. 6). In addition, the minimum peak HS date change is found during WW seasons (Table 3) because the snowpack would be scarce at those times, and there were no defined peaks (Fig. S4).

We determined the snow ablation sensitivity to temperature and precipitation change in response to different temperature increases at different elevations and during different compound season types. The results show there were small differences in the absolute snow ablation values in a warmer climate (Fig. 7). At low elevations, the average snow ablation sensitivity to temperature and precipitation change in all four compound seasons was 12 % ^∘C⁻¹ (Table 3). At mid and high elevations, the maximum snow ablation sensitivity to temperature and precipitation change was during dry seasons; WD seasons had a snow ablation sensitivity to temperature and precipitation change of 13 % ^∘C⁻¹ at mid elevations and 10 % ^∘C⁻¹ at high elevations. On the other hand, the minimum values for mid elevations were during WW seasons, when the snow ablation sensitivity to temperature and precipitation change was 8 % ^∘C⁻¹; the minimum values at high elevations were during CW seasons at 5 % ^∘C⁻¹.

Figure 6 Difference (days) from the baseline peak HS date at three different elevations and during the four different temperature (colors) and precipitation shifts (columns) for each season (boxes).

Table 3 Snow duration, snow ablation, and peak HS date sensitivity to temperature and precipitation change during the four different compound season types.

Figure 7 Absolute snow ablation values (centimeter per day) (y axis) at three different elevations during four different compound temperatures and precipitation for the baseline and different increments of temperature (x axis) seasons.

4.3 Spatial patterns

PCA reveals four Pyrenean sectors, namely northern–western (NW), northern–eastern (NE), southern–western (SW), and southern–eastern (SE). No remarkable differences between sectors are found in the relative importance of each compound season type in the snow sensitivity to temperature and precipitation change (Fig. 8). Snow sensitivity to temperature and precipitation change absolute values are generally lower at northern slopes (NW and NE) than at the southern slopes (SW and SE) (Figs. S7 and S8). In detail, seasonal HS ranged from −26 % ^∘C⁻¹ during CD seasons (NW) to −36 % ^∘C⁻¹ during WW (SE). Similarly, the maximum peak HS max sensitivity to temperature and precipitation change was at the SE sector during WW seasons (25 % ^∘C⁻¹) and the minimum was during CD seasons at the NW sector (15 % ^∘C⁻¹). The snow duration sensitivity to temperature and precipitation change increased during WW seasons, and the maximum changes were at the SE sector (−16 % ^∘C⁻¹); in contraposition, the lowest sensitivity to temperature and precipitation change is found at the NW sector during CD and CW seasons (−8 % ^∘C⁻¹ in both seasons). Snow ablation sensitivity to temperature and precipitation change increases towards the eastern Pyrenees, particularly during WD seasons (14 % ^∘C⁻¹ and 13 % ^∘C⁻¹ for the NE and SE sectors, respectively). Finally, no remarkable peak HS date differences are observed between sectors, and maximum values are found during CD and CW seasons when the peak HS date is anticipated to be ≥5 d per degree Celsius for all sectors.

Figure 8 Average snow sensitivity to temperature and precipitation change (y axis) grouped by sector (x axis), season type (color bars), and snow climate indicator (boxes).

5 Discussion

The spatial and temporal patterns of snow in the Pyrenees are highly variable, and climate projections indicate that extreme events will likely increase in future decades (Meng et al., 2022). Therefore, we analyzed factors that affect the snowpack sensitivity to temperature and precipitation change to gain insight into how future climate changes may affect the snow regime.

5.1 Snow sensitivity to temperature and precipitation change and its relationship with historical and future snow trends

5.1.1 Snow accumulation phase

The snow losses due to warming that we described here are mainly associated with increases in the rain : snowfall ratio (Fig. S9); changes in the snow onset and offset dates (Fig. S4); and increases in the energy available for snow ablation during the later months of the snow season, as has previously been reported in the literature (e.g., Pomeroy et al., 2015; Lynn et al., 2020; Jennings and Molotch, 2020). At high elevations, a trend of increasing precipitation (+10 %) could counterbalance temperature increases (<1 ^∘C; Fig. S5), consistent with the results previously reported for specific sites of the central Pyrenees (Izas, 2000 m; López-Moreno et al., 2008). Rasouli et al. (2014) also found that a 20 % increase in precipitation could compensate for a 2 ^∘C increase in temperature in subarctic Canada. A climate sensitivity analysis in the western Cascades (western USA) found that increases in precipitation due to warming modulated the snowpack accumulation losses by about 5 % 1^∘C⁻¹ (Minder, 2010). These results are consistent with recent data that examined snow above 1000 m in the Pyrenees, which found that an increase in the frequency of west circulation weather types since the 1980s increased the HS (Serrano-Notivoli et al., 2018; López-Moreno et al., 2020a), snow accumulation (Bonsoms et al., 2021a), and changes in winter snow days (Buisan et al., 2016). There are similar trends in the Alps, with an increase in extreme (exceeding the 100-year return level) snowfall above 3000 m during recent decades (Le Roux et al., 2021) and increases in extreme winter precipitation (Rajczak and Schär, 2017).

5.1.2 Snow ablation phase

Climate warming leads to a cascade of physical changes in the SEB that increase snow ablation near the 0 ^∘C isotherm. Overall, the snow ablation showed low to inexistent changes due to warming. A comparison between low and high elevations indicated slightly faster snow ablation at high elevations (Fig. 7). This higher rate of snow ablation per season at high elevations (which have deeper snowpacks) is probably because the snow there lasts until late spring, when more energy is available for snow ablation (Bonsoms et al., 2022). Temperature increase does not imply remarkable changes in snow ablation per season because warming decreases the magnitude of the snowpack (seasonal HS and peak HS max) and triggers an earlier onset of snowmelt (Wu et al., 2018). The earlier peak HS date (Table 3 and Fig. 6) implies lower rates of net shortwave radiation because snow melting starts earlier in warmer climates (Pomeroy et al., 2015), coinciding with the shorter days and lower solar zenith angle (Lundquist et al., 2013; Sanmiguel-Vallelado et al., 2022). Our results agree with the slow snow melt rates reported in the Northern Hemisphere from 1980 to 2017 (Wu et al., 2018). The results of previous studies were similar for subarctic Canada (Rasouli et al., 2014) and western USA snowpacks (Musselman et al., 2017b), but Arctic sites had faster melt rates (Krogh and Pomeroy, 2019).

5.1.3 Snow sensitivity to temperature and precipitation change and snowpack projections

Our results suggest that warming had a non-linear effect on snowpack reduction. Our largest snow losses were for seasonal HS when the temperature increased by 1 ^∘C above the baseline. At low and mid elevations, the average seasonal HS decrease was more than 40 % for all compound season types, and the maximum sensitivity was during WW seasons. Previous research in the Pyrenees and other mid-latitude mountain ranges reported similar results. A study in the central Pyrenees reported the peak SWE was 29 % ^∘C⁻¹, whereas snow season duration decreased by about 20 to 30 d at about 2000 m (López-Moreno et al., 2013a). The average peak HS max at high elevations in the Pyrenees (−16 % ^∘C⁻¹; Fig. 6 and Table 2) was similar to the average peak SWE sensitivity (−15 % ^∘C⁻¹) reported in the Iberian Peninsula mountains at 2500 m (Alonso-González et al., 2020b). These results are also consistent with climate projections for this mountain range. In particular, for a 2 ^∘C or more increase in temperature, the snow season declined by 38 % at the lowest ski resorts (∼1500 m) in the SE Pyrenees (Pons et al., 2015). However, high-emission climate scenarios projected an increase in the frequency and intensity of high snowfall at high elevations (López-Moreno et al., 2011a). Snow sensitivity in the easternmost areas could decline during the winter because of a trend of an increase of about 10 % in precipitation in this area (Amblar-Francés et al., 2020). Our projected changes in the Pyrenean snowpack dynamics are similar to the expected snow losses in other mountain ranges. For example, a study of the Atlas Mountains of northern Africa concluded that snowpack decreases were greater in the lowlands, and it projected seasonal SWE declines of 60 % under the RCP4.5 scenario and 80 % under the RCP8.5 scenario for the entire range (Tuel et al., 2022). A study in the Washington Cascades (western USA) found that snowpack decline was 19 % to 23 % per 1 ^∘C (Minder, 2010), similar to the values in the present study at high elevations. A study of the French Alps (Chartreuse, 1500 m) found that seasonal HS decreases of the order of 25 % for a 1.5 ^∘C increase and 32 % for a 2 ^∘C increase in global temperature above the pre-industrial years (Verfaille et al., 2018). A study of the Swiss Alps reported a snowpack decrease of about 15 % ^∘C⁻¹ (Beniston et al., 2003); in the same alpine country, another study predicted the seasonal HS will decrease by more than 70 % in massifs below 1000 m in all future climate projections (Marty et al., 2017). The largest snow reductions will likely occur during the periods between seasons (Steger et al., 2013; Marty et al., 2017). Nevertheless, at high elevations, snow climate projections found no significant trend for maximum HS until the end of the 21st century above 2500 m in the eastern Alps (Willibald et al., 2020), suggesting that internal climate variability is a major source of uncertainty of SWE projections at high elevations (Schirmer et al., 2022).

5.2 Influence of compound temperature and precipitation seasons

We found that the maximum sensitivities of seasonal HS and peak HS max to temperature and precipitation change were during WW seasons at low and mid elevations and during WD seasons at high elevations. Brown and Mote (2009) analyzed the sensitivity of snow to climate changes in the Northern Hemisphere and found maximal SWE sensitivities in mid-latitudinal maritime winter climate areas and minimal SWE sensitivities to temperature and precipitation change in dry and continental zones, consistent with our results. López-Moreno et al. (2017) also found greater decreases of SWE in wet and temperate Mediterranean ranges than in drier regions. Furthermore, Rasouli et al. (2022) studied the northern North American Cordillera and found higher snowpack sensitivities to temperature and precipitation change in wet basins than in dry basins. Our maximum snow ablation relative change over the baseline scenario occurred during WD seasons in accordance with Musselman et al. (2017b), who found a higher snow melt rate during dry years in the western USA. Low and mid elevations are highly sensitive to WW seasons because wet conditions favor decreases in the seasonal HS due to advection from sensible heat fluxes. The temperature in the Pyrenees is still cold enough to allow snowfall at high elevations during WW seasons, and for this reason we found maximal sensitivities to temperature and precipitation change during WD seasons. Reductions of snowfall in alpine regions can be compensated for in a warmer scenario because warm and wet snow is less susceptible to blowing wind transport and losses from sublimation (Pomeroy and Li, 2000; Pomeroy et al., 2015). During spring, snow runoff could also be greater in wet climates due to rain-on-snow events (Corripio and López-Moreno, 2017), coinciding with the availability of more energy for snow ablation.

5.3 Spatial and elevation factors controlling snow sensitivity to temperature and precipitation change

A comparison between Pyrenean sectors (Fig. 8) reveals no remarkable differences in the relative importance of each compound season type in the snow sensitivity to temperature and precipitation change. This is because by applying a joint-quantile approach for each massif and elevation, we are comparing similar climate seasons between sectors, regardless of the number of compound season types recorded in each massif during the baseline period (Figs. S1 and S3). The highest absolute snow sensitivity to temperature and precipitation change values is found in the SE Pyrenees. This is consistent with the snow accumulation and ablation patterns previously reported in this region (López-Moreno, 2005; Navarro-Serrano et al., 2017; Alonso-González et al., 2020b; Bonsoms et al., 2021a, b; Bonsoms et al., 2022). The Atlantic climate has a lower influence on the SE sector, and in situ observations indicated there was about half of the seasonal snow accumulation amounts as in northern and western areas at the same elevation (>2000 m; Bonsoms et al., 2021a). The snow in the SE Pyrenees is more sensitive to temperature and precipitation change because these massifs are exposed to higher turbulence and radiative heat fluxes (Bonsoms et al., 2022), and the 0 ^∘C isotherm is closer. Similar conclusions are found for low elevations where the results show an upward displacement of the snow line due to warming. Previous studies described the sensitivity of the snow pattern to elevation at specific stations of the central Pyrenees (López-Moreno et al., 2013b, 2017), Iberian Peninsula mountains (Alonso-González et al., 2020b), and other ranges such as the Cascades (Jefferson, 2011; Sproles et al., 2013), the Alps (Marty et al., 2017), and western USA (Pierce and Cayan, 2013; Musselman et al., 2017b). In these regions, the models suggest larger snowpack reductions due to warming at subalpine sites than at alpine sites (Jennings and Molotch, 2020) due to closer isothermal conditions (Brown and Mote, 2009; López-Moreno et al., 2017; Mote et al., 2018).

5.4 Environmental and socioeconomic implications

Our results indicated there will be an increase in snow ablation days and implied a disappearance of the typical sequence of snow accumulation and snow ablation seasons. Climate warming triggers the simultaneous occurrence of several periods of snowfall and melting, snow droughts during the winter, and ephemeral snowpacks between seasons. These expected decreases in snow will likely have important impacts on the ecosystem. During spring, a snow cover cools the soil (Luetschg et al., 2008), delays the initiation of freezing (Oliva et al., 2014), functions as a thick active layer (Hrbáček et al., 2016), and protects alpine rocks from exposure to solar radiation and high air temperatures (Magnin et al., 2017). Due to warming temperatures, the remaining glaciers in the Pyrenees are shrinking and are expected to disappear before the 2050s (Vidaller et al., 2021). The shallower snowpack that we identified in this work will increase the vulnerability of glaciers because snow has a higher albedo than dark ice and debris-covered glaciers and functions as a protective layer for glaciers (Fujita and Sakai, 2014).

The earlier onset of snowmelt suggested by our results, which is greater at low and mid elevations during WD seasons, is in line with previous global studies that have reported earlier streamflow due to earlier runoff dates (Adam et al., 2009; Stewart, 2009) and with a study of changes in the Iberian Peninsula River flows (Morán-Tejeda et al., 2014). Overall, our results are consistent with the slight decrease in the river peak flows that have occurred in the southern slopes of the Pyrenees since the 1980s (Sanmiguel-Vallelado et al., 2017). The reductions of seasonal HS that we identified suggest that snowmelt-dominated streamflows will likely shift to rainfall-dominated regimes. Although high-elevation meltwater might increase and contribute to earlier groundwater recharging (Evans et al., 2018), the increased evapotranspiration in the lowlands (Bonsoms et al., 2022) could counter this effect, so there is no net change in downstream areas (Stahl et al., 2010). Snow ephemerality triggers lower spring and summer flows (Barnett et al., 2005; Adam et al., 2009; Stahl et al., 2010) and has an impact on hydrological management strategies. Winter snow accumulation affects hydrological availability during the months when water and hydroelectric demands are higher. This is because reservoirs store water during periods of peak flows (winter and spring) and release water during the driest season in the lowlands (summer) (Morán-Tejeda et al., 2014). Recurrent snow-scarce seasons may intensify these hydrological impacts and lead to competition for water resources among different ecological and socioeconomic systems. The economic viability of mountain ski resorts in the Pyrenees depends on a regular snow cover (Gilaberte-Burdalo et al., 2014; Pons et al., 2015), but this is highly variable, especially at low and mid elevations. The expected increase in snow-scarce seasons that we identified here is consistent with climate projections for this region, which suggest that no Pyrenean ski resorts will be viable under the RCP8.5 scenario by the end of the 21st century (Spandre et al., 2019).

5.5 Limitations and uncertainties

The meteorological input data that we used to model snow were estimated for flat slopes, and the regionalization system we used was based on the SAFRAN system. According to this system, a mountain range is divided into massifs with homogeneous topography. The SAFRAN system has negative biases in shortwave radiation, a temperature precision of about 1 K, and biases in the accumulated monthly precipitation of about 20 kg m² (Vernay et al., 2022). Our estimates of snow sensitivity to temperature and precipitation change were based on the delta approach, which considers changes in temperature and precipitation based on climate projections for the Pyrenees (Amblar-Francés et al., 2020) but assumes that the meteorological patterns of the reference period will be constant over time. In this work we used a physical-based snow model, since it provides better results for future snow climate change estimations than degree-day models (Carletti et al., 2022). FSM2 is a physics-based model of intermediate complexity, and the estimates of snow densification are simpler than those from more complex models of snowpack. However, a more complex model does not necessarily perform better in terms of snowpack and runoff estimation (Magnusson et al., 2015). The FSM2 configuration implemented in this work includes snow meltwater retention, snowpack refreezing, and snow albedo based on snow age, which are the physical parameters included in the best-performing snow models according to Essery et al. (2013). Snow model sensitivity studies reveal that intermediate-complexity models exhibit similar snow depth accuracies as most complex multi-layer snow models, as well as robust performances across seasons (Terzago et al., 2020).

6 Conclusions

Our study assessed the impact of temperature and precipitation change on the Pyrenean snowpack during compound cold–hot and wet–dry seasons using a physical-based snow model that was forced by reanalysis data. We determined the snow sensitivity to temperature and precipitation change using five key indicators of snow accumulation and snow ablation. The lowest snow sensitivity to temperature and precipitation change was at high elevations of the NW Pyrenees and increased at lower elevations and in the SE slopes. An increase of 1 ^∘C at low- and mid-elevation regions led to remarkable decreases in the seasonal HS and snow duration. However, at high elevations precipitation plays a key role, and temperature is far from the isothermal 0 ^∘C during the middle of winter. In this region, a 10 % increase in precipitation over the eastern regions of this range, as suggested by the Spanish Meteorological Agency (AEMET), could compensate for temperature increases of the order of about < 1 ^∘C. The impact of climate warming depends on a combination of temperature and precipitation during compound seasons. Our analysis of seasonal HS and peak HS max indicated the greatest declines were during WW seasons, and the smallest declines were during CD seasons, independently of the Pyrenean sector. For snow duration, however, the highest (lowest) sensitivity to temperature and precipitation change is found during WD (CW) seasons. Similarly, snow ablation had slightly greater sensitivities to temperature and precipitation change during WD seasons in that snow ablation variation is less than 10 %, and the peak HS date occurred about 5 d earlier per degree Celsius. Our findings thus provide evidence that the Pyrenean snowpack is highly sensitive to climate warming and suggest that the snowpacks of other mid-latitude mountain ranges may also show a similar response to warming.