Future snowfall in the Alps: projections based on the EURO-CORDEX regional climate models

Abstract

Twenty-first century snowfall changes over the European Alps are assessed based on high-resolution regional climate model (RCM) data made available through the EURO-CORDEX initiative. Fourteen different combinations of global and regional climate models with a target resolution of 12 $\mathrm{km}$ and two different emission scenarios are considered. As raw snowfall amounts are not provided by all RCMs, a newly developed method to separate snowfall from total precipitation based on near-surface temperature conditions and accounting for subgrid-scale topographic variability is employed. The evaluation of the simulated snowfall amounts against an observation-based reference indicates the ability of RCMs to capture the main characteristics of the snowfall seasonal cycle and its elevation dependency but also reveals considerable positive biases especially at high elevations. These biases can partly be removed by the application of a dedicated RCM bias adjustment that separately considers temperature and precipitation biases.

Snowfall projections reveal a robust signal of decreasing snowfall amounts over most parts of the Alps for both emission scenarios. Domain and multi-model mean decreases in mean September–May snowfall by the end of the century amount to $-$25 and $-$45 % for representative concentration pathway (RCP) scenarios RCP4.5 and RCP8.5, respectively. Snowfall in low-lying areas in the Alpine forelands could be reduced by more than$-$80 %. These decreases are driven by the projected warming and are strongly connected to an important decrease in snowfall frequency and snowfall fraction and are also apparent for heavy snowfall events. In contrast, high-elevation regions could experience slight snowfall increases in midwinter for both emission scenarios despite the general decrease in the snowfall fraction. These increases in mean and heavy snowfall can be explained by a general increase in winter precipitation and by the fact that, with increasing temperatures, climatologically cold areas are shifted into a temperature interval which favours higher snowfall intensities. In general, percentage changes in snowfall indices are robust with respect to the RCM postprocessing strategy employed: similar results are obtained for raw, separated, and separated–bias-adjusted snowfall amounts. Absolute changes, however, can differ among these three methods.

Introduction

Snow is an important resource for the Alpine regions, be it for tourism, hydropower generation or water management (Abegg et al., 2007). According to the Swiss Federal Office of Energy hydropower generation accounts for approximately 55 % of the Swiss electricity production (SFOE, 2014). Consideration of changes in snow climatology needs to address aspects of both snow cover and snowfall. In the recent past, an important decrease in the mean snow cover depth and duration in the Alps was observed (e.g, Laternser and Schneebeli, 2003; Marty, 2008; Scherrer et al., 2004). Projections of future snow cover changes based on climate model simulations indicate a further substantial reduction (Schmucki et al., 2015a; Steger et al., 2013), which is strongly linked to the expected rise in temperatures (e.g. CH2011, 2011; Gobiet et al., 2014). On regional and local scales rising temperatures exert a direct influence on snow cover in two ways: first, total snowfall sums are expected to decrease as a result of a lower probability for precipitation to fall as snow, implying a decreasing overall snowfall fraction (ratio between solid and total precipitation). Second, snow on the ground is subject to faster and accelerated melt. These warming-induced trends might be modulated, for instance, by changes in atmospheric circulation patterns.

Although the snowfall fraction is expected to decrease during the 21st century (e.g. Räisänen, 2016), extraordinary snowfall events can still leave a trail of destruction. A recent example was the winter 2013–2014 with record-breaking heavy snowfall events along the southern rim of the European Alps (e.g. Techel et al., 2015). The catastrophic effects of heavy snowfall range from avalanches and floods to road or rail damage. In extreme cases these events can even result in the weight-driven collapse of buildings or loss of human life (Marty and Blanchet, 2011). Also, mean snowfall conditions, such as the mean number of snowfall days in a given period, can be of high relevance for road management (e.g. Zubler et al., 2015) or airport operation. Projections of future changes in snowfall, including mean and extreme conditions, are therefore highly relevant for long-term planning and adaptation purposes in order to assess and prevent related socio-economic impacts and costs.

Twenty-first century climate projections typically rely on comprehensive climate models. For large-scale projections, global climate models (GCMs) with a rather coarse spatial resolution of 100 $\mathrm{km}$ or more are used. To assess regional- to local-scale impacts, where typically a much higher spatial resolution is required, a GCM can be dynamically downscaled by nesting a regional climate model (RCM) over the specific domain of interest (Giorgi, 1990). In such a setup, the GCM provides the lateral and sea surface boundary conditions to the RCM. One advantage of climate models is their ability to estimate climate change in a physically based manner under different greenhouse gas (GHG) emission scenarios. With the Intergovernmental Panel on Climate Change's (IPCC) release of the Fifth Assessment Report (IPCC, 2013), the so-called representative concentration pathway (RCP) scenarios have been introduced (Moss et al., 2010). These specify GHG concentrations and corresponding emission pathways for several radiative forcing targets. To estimate inherent projection uncertainties, ensemble approaches employing different climate models, different greenhouse gas scenarios and/or different initial conditions are being used (e.g. Deser et al., 2012; Hawkins and Sutton, 2009; Rummukainen, 2010).

Within the last few years several studies targeting the future global and European snowfall evolution based on climate model ensembles were carried out (e.g. de Vries et al., 2013, 2014; Krasting et al., 2013; O'Gorman, 2014; Piazza et al., 2014; Räisänen, 2016; Soncini and Bocchiola, 2011). Most of these analyses are based on GCM output or older generations of RCM ensembles at comparatively low spatial resolution, which are not able to properly resolve snowfall events over regions with complex topography. New generations of high-resolution RCMs are a first step toward an improvement on this issue. This is in particular true for the most recent high-resolution regional climate change scenarios produced by the global CORDEX initiative (Giorgi et al., 2009) and its European branch EURO-CORDEX (Jacob et al., 2014). The present work aims to exploit this recently established multi-model archive with respect to future snowfall conditions over the area of the European Alps. It thereby complements the existing works of Piazza et al. (2014) and de Vries et al. (2014). These two works also exploit comparatively high-resolved RCM experiments but with a smaller focus domain in the case of Piazza et al. (2014; French Alps only) and based on a single-model ensemble with a comparatively small ensemble size (eight members) in the case of de Vries et al. (2014).

In general and on decadal to centennial timescales, two main drivers of future snowfall changes over the European Alps with competing effects on snowfall amounts are apparent from the available literature: (1) mean winter precipitation is expected to increase over most parts of the European Alps and in most EURO-CORDEX experiments (e.g. Rajczak et al., 2017; Smiatek et al., 2016), which in principle could lead to higher snowfall amounts. (2) Temperatures are projected to considerably rise throughout the year (e.g. Gobiet et al., 2014; Smiatek et al., 2016; Steger et al., 2013) with the general effect of a decreasing snowfall frequency and fraction, thus potentially leading to a reduction in overall snowfall amounts. Separating the above two competing factors is one of the targets of the current study. A potential complication is that changes in daily precipitation frequency (events with precipitation $>\mathrm{1}$ $\mathrm{mm}$ day${}^{-\mathrm{1}}$) and precipitation intensity (average amount on wet days) can change in a counteracting manner (e.g. Fischer et al., 2015; Rajczak et al., 2013) and that relative changes are not uniform across the event category (e.g.; Fischer and Knutti, 2016; Rajczak et al., 2017).

We here try to shed more light on these issues by addressing the following main objectives:

Snowfall separation on the RCM grid. Raw snowfall outputs are not available for all members of the EURO-CORDEX RCM ensemble. Therefore, an adequate snowfall separation technique, i.e. the derivation of snowfall amounts based on readily available daily near-surface air temperature and precipitation data, is required. Furthermore, we seek for a snowfall separation method that accounts for the topographic subgrid-scale variability in snowfall on the RCM grid.

Snowfall bias adjustment. Even the latest generation of RCMs is known to suffer from systematic model biases (e.g. Kotlarski et al., 2014). In GCM-driven setups as employed within the present work these might partly be inherited from the driving GCM. To remove such systematic biases in temperature and precipitation, a dedicated bias adjustment method is developed and employed in the present work. To assess its performance and applicability, different snowfall indices in the bias-adjusted and not-bias-adjusted output are compared against observation-based estimates.

Snowfall projections for the late 21st century. Climate change signals for various snowfall indices over the Alpine domain and for specific elevation intervals, derived by a comparison of 30-year control and scenario periods, are analysed under the assumption of the RCP8.5 emission scenario. In addition, we aim to identify and quantify the main drivers of future snowfall changes and, in order to assess emission scenario uncertainties, compare RCP8.5-based results with experiments assuming the more moderate RCP4.5 emission scenario. Snowfall projections are generally based on three different data sets: (1) raw RCM snowfall where available, (2) RCM snowfall separated from simulated temperature and precipitation, and (3) RCM snowfall separated from simulated temperature and precipitation and additionally bias-adjusted. While all three estimates are compared for the basic snowfall indices in order to assess the robustness of the projections, more detailed analyses are based on data set (3) only.

In addition and as preparatory analysis, we carry out a basic evaluation of RCM-simulated snowfall amounts. This evaluation, however, is subject to considerable uncertainties as a high-quality observation-based reference on the required spatial scale is not available, and the focus of the present work is on snowfall projections.

The article is structured as follows: Sect. 2 describes the data used and methods employed. In Sects. 3 and 4, results of the bias adjustment approach and snowfall projections for the late 21st century are shown, respectively. The latter are further discussed in Sect. 5 while overall conclusions and a brief outlook are provided in Sect. 6. Additional supporting figures are provided in the Supplement.

Data and methods

Observational data

To estimate the reference fine-scale snowfall, two gridded data sets, one for precipitation and one for temperature, derived from station observations and covering the area of Switzerland are used. Both data sets are available on a daily basis with a horizontal resolution of 2 $\mathrm{km}$ for the entire evaluation period (EVAL) 1971–2005 (see Sect. 2.3).

The gridded precipitation data set (RhiresD) represents a daily analysis based on a high-resolution rain-gauge network (MeteoSwiss, 2013a) consisting of more than 400 stations that have a balanced distribution in the horizontal but under-represent high altitudes (Frei and Schär, 1998; Isotta et al., 2014). Despite the data set's resolution of 2 $\mathrm{km}$, the effective grid resolution as represented by the mean inter-station distance is about 15–20 $\mathrm{km}$ and thus comparable to the nominal resolution of the available climate model data (see Sect. 2.2). The data set has not been corrected for the systematic measurement bias of rain gauges (e.g. Neff, 1977; Sevruk, 1985; Yang et al., 1999).

The gridded near-surface air temperature (from now on simply referred to as temperature) data set (TabsD) utilises a set of approximately 90 homogeneous long-term station series (MeteoSwiss, 2013b). Despite the high quality of the underlying station series, errors might be introduced by unresolved scales, an uneven spatial distribution and interpolation uncertainty (Frei, 2014). The unresolved effects of land cover or local topography, for instance, probably lead to an underestimation of spatial variability. Also note that, while RhiresD provides daily precipitation sums aggregated from 06:00 to 06:00 UTC of the following day, TabsD is a true daily temperature average from midnight UTC to midnight UTC. Due to a high temporal autocorrelation of daily mean temperature, this slight inconsistency in the reference interval of the daily temperature and precipitation grids is expected to not systematically influence our analysis.

Overview on the full and the reduced EURO-CORDEX simulations employed in this study. All experiments were realised on a grid covering the European domain with a horizontal resolution of approximately 12 $\mathrm{km}$ (EUR-11) and were run for the emission scenarios RCP4.5 and RCP8.5. A subset of nine simulations provides raw snowfall, i.e. snowfall flux in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ s, as an output variable (marked by an “X”). The reduced model set excludes experiments that are subject to serious shortcomings (see Supplement, Part B). For full institutional names the reader is referred to the official EURO-CORDEX website www.euro-cordex.net. Note that the GCM name (column “GCM”) includes the name of the institution that carried out the specific GCM run. Note that the EC-EARTH-driven experiments partly employ different realisations of the GCM run; i.e. they explicitly sample the influence of internal climate variability in addition to model uncertainty.

RCM	Institute	GCM	Abbreviation	Raw snowfall	Part of reduced
			used in this work		model set
ALADIN53	CNRM	CNRM-CERFACS-CNRM-CM5	CNRM–ALADIN	–	X
CCLM4-8-17	CLMcom/BTU	CNRM-CERFACS-CNRM-CM5	CNRM–CCLM	–	X
CCLM4-8-17	CLMcom/BTU	ICHEC-EC-EARTH${}^{\mathrm{d}}$	EC-EARTH–CCLM	–	X
CCLM4-8-17	CLMcom/ETH	MOHC-HadGEM2-ES	HadGEM2–CCLM	–	X
CCLM4-8-17	CLMcom/BTU	MPI-M-MPI-ESM-LR	MPI-ESM–CCLM	–	X
HIRHAM5	DMI	ICHEC-EC-EARTH${}^{\mathrm{c}}$	EC-EARTH–HIRHAM	X	X
RCA4	SMHI	CNRM-CERFACS-CNRM-CM5	CNRM–RCA	X	X
RCA4	SMHI	ICHEC-EC-EARTH${}^{\mathrm{d}}$	EC-EARTH–RCA	X	X
RCA4	SMHI	MOHC-HadGEM2-ES	HadGEM2–RCA	X	X
RCA4	SMHI	IPSL-IPSL-CM5A-MR	IPSL–RCA	X	X
RCA4	SMHI	MPI-M-MPI-ESM-LR	MPI-ESM–RCA	X	X
REMO2009	MPI-CSC	MPI-M-MPI-ESM-LR${}^{\mathrm{a}}$	MPI-ESM–REMO	X	X
RACMO22E	KNMI	ICHEC-EC-EARTH${}^{\mathrm{b}}$	EC-EARTH–RACMO	X	–
WRF331F	IPSL-INERIS	IPSL-IPSL-CM5A-MR	IPSL–WRF	X	–

${}^{\mathrm{a}}$ r1i1p1 realisation. ${}^{\mathrm{b}}$ r1i1p1 realisation. ${}^{\mathrm{c}}$ r3i1p1 realisation. ${}^{\mathrm{d}}$ r12i1p1 realisation.

In addition to the gridded temperature and precipitation data sets and in order to validate simulated raw snowfall amounts, station-based observations of fresh snow sums (snow depth) at daily resolution from 29 stations in Switzerland with data available for at least 80 % of the evaluation period 1971–2005 are employed.

GTOPO30 topography (https://lta.cr.usgs.gov/GTOPO30) aggregated to the EUR-11 (0.11${}^{\circ }$) RCM grid. The coloured area shows the Alpine domain used for the assessment of snowfall projections. The bold black outline marks the Swiss subdomain used for the assessment of the bias adjustment approach.

[Figure omitted. See PDF]

Climate model data

In terms of climate model data we exploit a recent ensemble of regional climate projections made available by EURO-CORDEX (www.euro-cordex.net), the European branch of the World Climate Research Programme's CORDEX initiative (www.cordex.org; Giorgi et al., 2009). RCM simulations for the European domain were run at a resolution of approximately 50 $\mathrm{km}$ (EUR-44) and 12 $\mathrm{km}$ (EUR-11) with both reanalysis boundary forcing (Kotlarski et al., 2014; Vautard et al., 2013) and GCM forcing (Jacob et al., 2014). We here disregard the reanalysis-driven experiments and employ the GCM-driven simulations only. These include historical control simulations and future projections based on RCP greenhouse gas and aerosol emission scenarios. We employ daily averaged model output of GCM-driven EUR-11 simulations that were available in December 2016 and for which control, RCP4.5 and RCP8.5 runs are provided. We hence exclusively focus on the higher-resolved EUR-11 simulations and disregard the coarser EUR-44 ensemble due to the apparent added value of the EUR-11 ensemble with respect to regional-scale climate features in the complex topographic setting of the European Alps (e.g. Giorgi et al., 2016; Torma et al., 2015). Out of the entire set of available EURO-CORDEX simulations, several GCM–RCM chains were either completely or partially removed from our analysis, resulting in a full set of 14 GCM–RCM combinations and a reduced set of 12 combinations only (Table 1). Reasons for removal are the existence of several realisations (MPI-ESM–REMO; only one realisation has been used) and simulation deficiencies that potentially affect our analysis (HadGEM2–RACMO, EC-EARTH–RACMO and IPSL–WRF). Further details on these issues are provided in the Supplement, Part B.

Analysed snowfall indices. The last column indicates the threshold value in the CTRL period for considering a grid cell in the climate changes analysis (grid cells with smaller values are skipped for the respective analysis): the first number is the threshold for monthly analyses and the second number is the threshold for seasonal analysis.

Index name	Abbreviation	Unit	Definition	Threshold for monthly/
				seasonal analysis
Mean snowfall	$S_{\text{mean}}$	mm	(Spatio-)temporal mean snowfall in mm snow water
			equivalent (only “mm” hereafter).	1 $\mathrm{mm}$/10 $\mathrm{mm}$
Heavy snowfall	$S_{\mathrm{q}\mathrm{99}}$	$\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$	Grid-point-based 99 % all-day-snowfall percentile.	1 $\mathrm{mm}$/1 $\mathrm{mm}$
Max. 1-day snowfall	$S_{\mathrm{1}\mathrm{d}}$	$\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$	Mean of each season's or month's maximum 1-day snowfall.	1 $\mathrm{mm}$/1 $\mathrm{mm}$
Snowfall frequency	$S_{\text{freq}}$	%	Percentage of days with snowfall $S>\mathrm{1}$ $\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$
			within a specific time period.	1 %/1 %
Snowfall intensity	$S_{\text{int}}$	$\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$	Mean snowfall intensity at days with snowfall $S>\mathrm{1}$ $\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$
			within a specific time period.	S${}_{\text{freq}}$ threshold passed
Snowfall fraction	$S_{\text{frac}}$	%	Percentage of total snowfall, $S_{\text{tot}}$, on total precipitation, $P_{\text{tot}}$,
			within a specific time period.	1 %/1 %

It is important to note that each of the RCMs considered uses an individual topography field. Model topographies for a given grid cell might therefore considerably differ from each other and also from the observation-based orography. It is therefore not meaningful to compare snowfall values at individual grid cells since the latter might be situated at different elevations. For this reason most analyses of the present work were carried out as a function of elevation, i.e. by averaging quantities over distinct elevation intervals.

Analysis domain and periods

The arc-shaped European Alps – with a West–East extent of roughly 1200 $\mathrm{km}$, a total of area 190 000 ${\mathrm{km}}^{\mathrm{2}}$ and a peak elevation of 4810 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Mont Blanc) – are the highest and most prominent mountain range which is entirely situated in Europe. In the present work, two different analysis domains are used. The evaluation of the bias adjustment approach depends on the observational data sets RhiresD and TabsD (see Sect. 2.1). As these cover Switzerland only, the evaluation part of the study (Sect. 3) is constrained to the Swiss domain (Fig. 1, bold line). For the analysis of projected changes in different snowfall indices (Sects. 4 and 5) a larger domain covering the entire Alpine crest with its forelands is considered (Fig. 1, coloured region).

Our analysis is based on three different time intervals. The evaluation period 1971–2005 is used for the calibration and validation of the bias adjustment approach. Future changes in snowfall indices are computed by comparing a present-day control period (1981–2010, CTRL) to a future scenario period at the end of the 21st century (2070–2099, SCEN). For all periods (EVAL, CTRL and SCEN), the summer months June, July and August (JJA) are excluded from the analysis. In addition to seasonal mean snowfall conditions, i.e. averages over the 9-month period from September to May, we also analyse the seasonal cycle of individual snowfall indices at monthly resolution.

Analysed snowfall indices and change signals

A set of six different snowfall indices is considered (Table 2). Mean snowfall ($S_{\text{mean}}$) refers to the (spatio-)temporally averaged snowfall amount in mm snow water equivalent (SWE; note that from this point on we will use the unit “mm” as a synonym for “mm SWE” as the unit of several snowfall indices). The two indices heavy snowfall ($S_{\mathrm{q}\mathrm{99}}$) and maximum 1-day snowfall (S${}_{\mathrm{1}\mathrm{d}}$) allow the assessment of projected changes in heavy snowfall events and amounts. $S_{\mathrm{1}\mathrm{d}}$ is derived by averaging maximum 1-day snowfall amounts over all individual months/seasons of a given time period (i.e. by averaging 30 maximum values in the case of the CTRL and SCEN period), while $S_{\mathrm{q}\mathrm{99}}$ is calculated from the grid-point-based 99th all-day-snowfall percentile of the daily probability density function (PDF) for the entire time period considered. We use all-day percentiles as the use of wet-day percentiles leads to conditional statements that are often misleading (see the analysis in Schär et al., 2016). Note that the underlying number of days differs for seasonal (September–May) and monthly analyses. Snowfall frequency ($S_{\text{freq}}$) and mean snowfall intensity ($S_{\text{int}}$) are based on a snow-day threshold of 1 $\mathrm{mm}$ day${}^{-\mathrm{1}}$ and provide additional information about the distribution and magnitude of snowfall events, while the snowfall fraction ($S_{\text{frac}}$) describes the ratio of solid precipitation to total precipitation. As climate models tend to suffer from too-high occurrence of drizzle and as small precipitation amounts are difficult to measure, daily precipitation values smaller than or equal to 0.1 $\mathrm{mm}$ were set to zero in both the observations and the simulations prior to the remaining analyses.

Snowfall ratios for the binary and Richards snow fractionation method. Ratios represent the quotient of the snowfall as estimated by the respective method and the snowfall as estimated by the subgrid method. Ratios are valid at the coarse-resolution grid (12 $\mathrm{km}$). (a) Ratios for mean snowfall, $S_{\text{mean}}$. (b) Ratios for heavy snowfall, $S_{\mathrm{q}\mathrm{99}}$. Ratio means were derived after averaging the corresponding snowfall index for 250 m elevation intervals in Switzerland while the ratio spread represents the minimum and maximum grid-point-based ratios in the corresponding elevation interval. This analysis is entirely based on the observational data sets TabsD and RhiresD.

[Figure omitted. See PDF]

Projections are assessed by calculating two different types of changes between the CTRL and the SCEN period. The absolute change signal ($\mathrm{\Delta }$) of a particular snowfall index $X$ (see Table 2) can be calculated as $${\displaystyle }{\displaystyle }\mathrm{\Delta }X=X_{\text{SCEN}}-X_{\text{CTRL}}$$ and the relative change signal ($\mathit{\delta }$) which describes the change in the snowfall index as a percentage of its CTRL period value: $${\displaystyle }{\displaystyle }\mathit{\delta }X=\left({\displaystyle \frac{X_{\text{SCEN}}}{X_{\text{CTRL}}}},-,\mathrm{1}\right)\cdot \mathrm{100}.$$ To prevent erroneous data interpretation due to possibly large relative changes in small CTRL values, certain grid boxes were masked out before calculating and averaging the signal of change. This filtering was done by setting threshold values for individual indices and statistics (Table 2).

Separating snowfall from total precipitation

Due to (a) the lack of a gridded observational snowfall data set and (b) the fact that not all EURO-CORDEX RCMs provide raw snowfall as an output variable, a method to separate solid precipitation from total precipitation depending on near-surface temperature conditions is developed. The simplest approach to separate snowfall from total precipitation is to fractionate the two phases in a binary manner by applying a constant snow fractionation temperature (e.g. de Vries et al., 2014; Schmucki et al., 2015a; Zubler et al., 2014). More sophisticated methods estimate the snow fraction $f_{\mathrm{s}}$ dependence on air temperature with linear or logistic relations (e.g. Kienzle, 2008; McAfee et al., 2014). In our case, the different horizontal resolutions of the observational (high resolution of 2 $\mathrm{km}$) and simulated (coarser resolution of 12 $\mathrm{km}$) data sets further complicate a proper comparison of the respective snowfall amounts. Thus, we explicitly analysed the snowfall amount dependency on the grid resolution and exploited possibilities for including subgrid-scale variability in snowfall separation. This approach is important as especially in Alpine terrain a strong subgrid-scale variability in near-surface temperatures due to orographic variability has to be expected, with corresponding effects on the subgrid-scale snowfall fraction.

For this preparatory analysis, which is entirely based on observational data, a reference snowfall is derived. It is based on the approximation of snowfall by application of a fixed temperature threshold to daily total precipitation amounts on the high-resolution observational grid (2 $\mathrm{km}$) and will be termed the subgrid method hereafter: first, the daily snowfall $S^{\prime }$ at each grid point of the observational data set at high resolution (2 $\mathrm{km}$) is derived by applying a snow fractionation temperature $T^{\ast }=\mathrm{2}$ ${}^{\circ }$C. The whole daily precipitation amount $P^{\prime }$ is accounted for as snow $S^{\prime }$ (i.e. $f_{\mathrm{s}}=\mathrm{100}$ %) for days with daily mean temperature $T^{\prime }\le T^{\ast }$. For days with $T^{\prime }>T^{\ast }$, $S^{\prime }$ is set to zero and $P^{\prime }$ is attributed as rain (i.e. $f_{\mathrm{s}}=\mathrm{0}$ %). This threshold approach with a fractionation temperature of 2 ${}^{\circ }$C corresponds to the one applied in previous works, and results appear to be in good agreement with station-based snowfall measurements (e.g. Zubler et al., 2014). The coarse grid (12 $\mathrm{km}$) reference snowfall $S_{\text{SG}}$ is determined by averaging the sum of separated daily high-resolution $S^{\prime }$ over all $n$ high-resolution grid points $i$ located within a specific coarse grid point $k$. That is, at each coarse grid point $k$, $${\displaystyle }{\displaystyle }{S}_{\text{SG}}={\displaystyle \frac{\mathrm{1}}{n}}\cdot {\sum }_{i=\mathrm{1}}^n{P}_i^{\prime }\left[T_i^{\prime },\le ,T^{\ast }\right]={\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n{S}_i^{\prime }.$$ For comparison, the same binary fractionation method with a temperature threshold of $T^{\ast }=\mathrm{2}$ ${}^{\circ }$C is directly applied to the coarse 12 $\mathrm{km}$ grid (binary method). For this purpose, total precipitation $P^{\prime }$ and daily mean temperature $T^{\prime }$ of the high-resolution data are conservatively remapped to the coarse grid leading to $P$ and $T$, respectively. Compared to the subgrid method, the binary method neglects any subgrid-scale variability in the snowfall fraction. As a result, the binary method underestimates $S_{\text{mean}}$ and overestimates $S_{\mathrm{q}\mathrm{99}}$ for most elevation intervals (Fig. 2). The underestimation of $S_{\text{mean}}$ can be explained by the fact that even for a coarse grid temperature above $T^{\ast }$, individual high-elevation subgrid cells (at which $T^{\prime }\le T^{\ast }$) can receive substantial snowfall amounts. As positive precipitation–elevation gradients can be assumed for most parts of the domain (larger total precipitation at high elevations; see e.g. Kotlarski et al., 2012 and 2015, for an Alpine-scale assessment), the neglect of subgrid-scale snowfall variation in the binary method hence leads to a systematic underestimation of mean snowfall compared to the subgrid method. Furthermore, following O'Gorman (2014), heavy snowfall events are expected to occur in a narrow temperature range below the rain–snow transition. As the binary method in these temperature ranges always leads to a snowfall fraction of 100 %, too-large $S_{\mathrm{q}\mathrm{99}}$ values would result.

To take into account these subgrid-scale effects, a more sophisticated approach – referred to as the Richards method – is developed here. This method is based upon a generalised logistic regression (Richards, 1959). Here, we apply this regression to relate the surface temperature $T$ to the snow fraction $f_{\mathrm{s}}$ by accounting for the topographic subgrid-scale variability. At each coarse grid point $k$, the Richards-method-based snowfall fraction $f_{\text{s, RI}}$ for a given day is hence computed as follows: $${\displaystyle }{f}_{\text{s,RI}}\left(T_k\right)={\displaystyle \frac{\mathrm{1}}{{\left[\mathrm{1},+,C_k,\cdot ,e^{D_k\cdot (T_k-T^{\ast })}\right]}^{\frac{\mathrm{1}}{C_k}}}},$$ with $C$ as the point of inflexion (denoting the point with the largest slope) and $D$ as the growth rate (reflecting the mean slope). $T_k$ is the daily mean temperature of the corresponding coarse grid box $k$ and $T^{\ast }=\mathrm{2}$ ${}^{\circ }$C is the snow fractionation temperature. First, we estimate the two parameters $C$ and $D$ of Eq. (4) for each single coarse grid point $k$ by minimising the least-squares distance to the $f_{\mathrm{s}}$ values derived by the subgrid method via the reference snowfall $S_{\text{SG}}$ (local fit). Second, $C$ and $D$ are expressed as a function of the topographic SD ${\mathit{\sigma }}_h$ of the corresponding coarse-resolution grid point only (Fig. S1 in the Supplement; global fit). This makes it possible to define empirical functions for both $C$ and $D$ that can be used for all grid points $k$ in the Alpine domain and that depend on ${\mathit{\sigma }}_h$ only. $$\begin{array}{}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{h,k}=\sqrt{{\displaystyle \frac{{\sum }_i^n{\left(h_i,-,\stackrel{\mathrm{‾}}{h_k}\right)}^{\mathrm{2}}}{n-\mathrm{1}}}}& {\displaystyle }& {\displaystyle }{C}_k={\displaystyle \frac{\mathrm{1}}{(E-{\mathit{\sigma }}_{h,k}\cdot F)}}& {\displaystyle }& {\displaystyle }{D}_k=G\cdot {\mathit{\sigma }}_{h,k}^{-H}\end{array}$$ Through a minimisation of the least-squares differences, the constant parameters in Eqs. (6) and (7) are calibrated over the domain of Switzerland and using daily data from the period from September to May 1971–2005, leading to values of $E=\mathrm{1.148336}$, $F=\mathrm{0.000966}$ ${\mathrm{m}}^{-\mathrm{1}}$, $G=\mathrm{143.84113}$ ${}^{\circ }$C${}^{-\mathrm{1}}$ and $H=\mathrm{0.8769335}$. Note that ${\mathit{\sigma }}_h$ is sensitive to the resolution of the two grids to be compared (see Eq. 5). It is a measure for the uniformity of the underlying topography and has been computed based on the high-resolution GTOPO30 digital elevation model (https://lta.cr.usgs.gov/GTOPO30) aggregated to a regular grid of 1.25 $\mathrm{arcs}$ (about 2 $\mathrm{km}$) which reflects the spatial resolution of the observed temperature and precipitation grids (see Sect. 2.1). Small values of ${\mathit{\sigma }}_h$ indicate a low subgrid-scale topographic variability, such as in the Swiss lowlands, while high values result from non-uniform elevation distributions, such as in areas of inner Alpine valleys. ${\mathit{\sigma }}_h$ as derived from GTOPO30 might be different from the subgrid-scale topographic variance employed by the climate models themselves, which is not relevant here as only grid-cell-averaged model output is analysed while ${\mathit{\sigma }}_h$ is regarded as a proper estimate of subgrid-scale variability.

Figure S1c provides an example of the relation between daily mean temperature $T$ and daily snow fraction $f_{\mathrm{s}}$ for grid cells with topographical SDs of 50 and 500 $\mathrm{m}$, respectively. The snowfall amount $S_{\text{RI}}$ for a particular day and a particular coarse grid box is finally obtained by multiplying the corresponding $f_{\text{s, RI}}$ and $P$ values. A comparison with the subgrid method yields very similar results. For both indices $S_{\text{mean}}$ and $S_{\mathrm{q}\mathrm{99}}$, mean ratios across all elevation intervals are close to 1 (Fig. 2). At single grid points, maximum deviations are not larger than $\mathrm{1}\pm \mathrm{0.1}$. Note that for this comparison calibration and validation period are identical (EVAL period). Based on this analysis, it has been decided to separate snowfall according to the Richards method throughout this work in both the observations and in the RCMs. The observation-based snowfall estimate obtained by applying the Richards method to the observational temperature and precipitation grids after spatial aggregation to the 0.11${}^{\circ }$ RCM resolution will serve as reference for the RCM bias adjustment and will be referred to as the “reference” hereafter. One needs to bear in mind that the parameters $C$ and $D$ of the Richards method were fitted for the Swiss domain only and were later applied to the entire Alpine domain (see Fig. 1).

Bias adjustment approach

Previous work has revealed partly substantial temperature and precipitation biases of the EURO-CORDEX RCMs over the Alps (e.g. Kotlarski et al., 2014; Smiatek et al., 2016), and one has to expect that the separated snowfall amounts are biased too. This would especially hamper the interpretation of absolute climate change signals of the considered snow indices. We therefore explore possibilities to bias-adjust the simulated snowfall amounts and to directly integrate this bias adjustment into the snowfall separation framework of Sect. 2.5. Note that we deliberately employ the term bias adjustment as opposed to bias correction to make clear that only certain aspects of the snowfall climate are adjusted and that the resulting data set might be subject to remaining inaccuracies.

A simple two-step approach that separately accounts for precipitation and temperature biases and their respective influence on snowfall is chosen. The separate consideration of temperature and precipitation biases allows for a more physically based bias adjustment of snowfall amounts: due to the temperature dependency of snowfall occurrence, snowfall biases of a given climate model cannot be expected to remain constant under changing climate conditions. For instance, a climate model with a given temperature bias might pass the snow–rain temperature threshold earlier or later than in reality during the general warming process. Hence, traditional bias adjustment approaches based only on a comparison of observed and simulated snowfall amounts in the historical climate would possibly fail due to a non-stationary bias structure. The bias adjustment is calibrated in the EVAL period for each individual GCM–RCM chain and over the region of Switzerland, and it is then applied to both the CTRL and SCEN period of each chain and for the entire Alpine domain. To be consistent in terms of horizontal grid spacing, the observational data sets RhiresD and TabsD (see Sect. 2.1) are conservatively regridded to the RCM resolution beforehand.

In a first step, total simulated precipitation is adjusted by introducing an elevation-dependent adjustment factor which adjusts precipitation biases regardless of temperature. For this purpose, mean precipitation ratios (RCM simulation divided by observational analysis) for 250 $\mathrm{m}$ elevation intervals are calculated (Fig. S3). An almost linear relationship of these ratios with elevation is found. Thus, a linear regression between the intervals from 250 to 2750 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ is used for each model chain separately to estimate a robust adjustment factor. As the total Alpine area above 2750 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ is small (Sect. 2.1; Fig. S2) and biases are subject to a considerable sampling uncertainty, these elevations are not considered in the regression. Overall the fits are surprisingly precise except for the altitude bins above 2000 $\mathrm{m}$ (Fig. S3). The precipitation adjustment factors $P_{\text{AF}}$ for a given elevation are then obtained as the inverse of the fitted precipitation ratios. Multiplying simulated precipitation $P$ with $P_{\text{AF}}$ for the respective model chain and elevation results in the adjusted precipitation: $${\displaystyle }{\displaystyle }{P}_{\text{adj}}=P\cdot P_{\text{AF}}.$$ For a given GCM–RCM chain and for each elevation interval, the spatially and temporally averaged corrected total precipitation $P_{\text{adj}}$ approximately corresponds to the observation-based estimate in the EVAL period.

In the second step of the bias adjustment procedure, temperature biases are accounted for. For this purpose the initial snow fractionation temperature $T^{\ast }=\mathrm{2}$ ${}^{\circ }$C of the Richards separation method (see Sec 2.5) is shifted to the value $T_{\mathrm{a}}^{\ast }$ for which the spatially (Swiss domain) and temporally (September to May) averaged simulated snowfall amounts for elevations below 2750 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ match the respective observation-based reference (see above). Compared to the adjustment of total precipitation, $T_{\mathrm{a}}^{\ast }$ is chosen independent of elevation but separately for each GCM–RCM chain in order to avoid over-parameterisation and to not over-interpret the elevation dependency of mean snowfall in the snowfall reference grid. After this second step of the bias adjustment, the spatially and temporally averaged simulated snowfall amounts below 2750 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ match the reference by definition. Hence, the employed simple bias adjustment procedure adjusts domain-mean snowfall biases averaged over the entire season from September to May. It does not, however, correct for biases in the spatial snowfall pattern, in the seasonal cycle or in the temporal distribution of daily values. Note that, as the underlying high-resolution data sets are available over Switzerland only, the calibration of the bias adjustment methodology is correspondingly restricted, but the adjustment is then applied to the whole Alpine domain. This approach is justified as elevation-dependent mean winter precipitation and temperature biases of the RCMs employed – assessed by comparison against the coarser-resolved EOBS reference data set (Haylock et al., 2008) – are very similar for Switzerland and for the entire Alpine analysis domain (Figs. S4 and S5).

Evaluation

RCM raw snowfall

We first carry out an illustrative comparison of RCM raw snowfall amounts (only for those simulations that directly provide snowfall flux) against station observations of snowfall in order to determine whether the simulated RCM snowfall climate contains valid information despite systematic biases. To this end, simulated raw snowfall amounts of seven EURO-CORDEX simulations of the reduced set (Table 1) averaged over 250 $\mathrm{m}$ elevation intervals and over the range 950–1650 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ are compared against observations of measured fresh snow sums from 29 MeteoSwiss stations (Sect. 2.1). For this purpose a mean snow density of 100 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ for the conversion from measured snow depth to water equivalent is assumed. Note that this simple validation is subject to considerable uncertainties as it does not explicitly correct for the scale and elevation gap between grid-cell-based RCM output and single-site observations. Especially in complex terrain and for exposed sites, point measurements of snow depth might be unrepresentative for larger-scale conditions (e.g. Grünewald and Lehning, 2015). Also, the conversion from snow depth to snow water equivalent is of approximate nature only, and fresh snow sums might furthermore misrepresent true snowfall in cases where snowmelt or snow drift occurs between two snow depth readings.

At low elevations simulated mean September–May raw snowfall sums match the observations well while differences are larger aloft (Fig. 3a). The positive bias at high elevations might arise from the fact that (the very few) observations were made at specific locations while simulated grid-point values of the corresponding elevation interval might be located in different areas of Switzerland. It might also be explained by positive RCM precipitation and negative RCM temperature biases at high elevations of the Alps (e.g. Kotlarski et al., 2015). Also note that, in general, the total high-elevation area of the Alpine analysis domain is small and elevations above 2500 $\mathrm{m}$ represent less than 5 % of the total area (Fig. S2). Both model-based and observation-based estimates for high elevations are hence subject to a considerable sampling uncertainty and are likely to be less robust than estimates for lower elevations.

Comparison of measured fresh snow sums of 29 MeteoSwiss stations (red) against simulated RCM raw snowfall in Switzerland (green) and against the 2 $\mathrm{km}$ reference snowfall grid obtained by employing the subgrid method (black) in the EVAL period 1971–2005. (a) Mean September–May snowfall vs. elevation. Both the simulation data (green) and the reference data (black) are based on the spatio-temporal mean of 250 $\mathrm{m}$ elevation ranges and plotted at the mean elevation of the corresponding interval. (b) Seasonal September–May snowfall cycle for the elevation interval 950 to 1650 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ Simulated multi-model means and spreads are based on a subset of seven EURO-CORDEX simulations of the reduced model set providing raw snowfall as an output variable (see Table 1).

[Figure omitted. See PDF]

At lower elevations, the station network is geographically more balanced and the observations are probably more representative of the respective elevation interval. Despite a clear positive snowfall bias in midwinter, the RCMs are generally able to reproduce the mean seasonal cycle of snowfall for elevations between 950 and 1650 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 3b). The fact that the major patterns of both the snowfall–elevation relationship and the mean seasonal snowfall cycle are well represented indicates the general and physically consistent applicability of RCM output to assess future changes in mean and heavy Alpine snowfall. However, substantial biases in snowfall amounts are apparent and a bias adjustment of simulated snowfall seems to be required prior to the analysis of climate change signals of individual snowfall indices.

Evaluation of the reference snowfall

The snowfall separation employing the Richards method (Sect. 2.5) and, as a consequence, also the bias adjustment (Sect. 2.6) make use of the 2 $\mathrm{km}$ reference snowfall grid derived by employing the subgrid method on the observed temperature and precipitation grids. Hence, the final results of this study could to some extent be influenced by inaccuracies and uncertainties of the reference snowfall grid itself. In order to assess the quality of the latter and in the absence of a further observation-based reference, we here present an approximate evaluation.

First, the reference snowfall grid is evaluated against fresh snow sums at the 29 Swiss stations that are also used for evaluating RCM raw snowfall. Note the limitations of such a comparison as outlined in Sect. 3.1. The comparison of black and red markers and lines in Fig. 3 indicates a good agreement of mean snowfall at individual elevation intervals (left panel) as well as for the mean annual cycle of snowfall at medium elevations (right panel). The reference snowfall grid is obviously a good approximation of site-scale fresh snow sums. Note that similarly to the RCM raw snowfall evaluation, all 2 $\mathrm{km}$ reference snowfall grid cells in the respective elevation interval are considered. The good agreement, however, still holds if only those 2 $\mathrm{km}$ grid cells covering the 29 site locations are considered (not shown here).

Overview of bias adjustment. (a) Elevation-dependent total precipitation adjustment factors, $P_{\text{AF}}$, for the full model set (14 GCM–RCM chains). (b) Scatter plot of mean September to May temperature biases (RCM simulation minus observational analysis) vs. adjusted snow fractionation temperatures, $T_{\mathrm{a}}^{\ast }$.

[Figure omitted. See PDF]

Second, both the 2 $\mathrm{km}$ reference snowfall grid and the 0.11${}^{\circ }$ reference snowfall grid obtained by employing the Richards method to aggregated temperature and precipitation values (Sect. 2.5) are compared against the gridded HISTALP data set of solid precipitation (Chimani et al., 2011). The latter is provided at a monthly resolution on a 5 min $\times$ 5 min grid covering the Greater Alpine Region. It is based on monthly snowfall fraction estimates that are used to scale a gridded data set of total precipitation. The comparison of the three data sets for the region of Switzerland (for which the 2 $\mathrm{km}$ reference snowfall is available) in the EVAL period 1971–2005 yields an approximate agreement of both the magnitude of mean winter snowfall and its spatial pattern (Fig. S6). The three data sets differ with respect to their spatial resolution but all show a clear dependency of snowfall on topography and mean September–May snowfall sums above 1000 $\mathrm{mm}$ over most parts of the Alpine ridge. Climatologically warm and dry valleys, on the other hand, are represented by minor snowfall amounts of less than 400 $\mathrm{mm}$ only.

As mentioned before these evaluations of the reference snowfall grid are subject to uncertainties and, furthermore, only cover mean snowfall amounts. However, they provide basic confidence in the applicability of the reference snowfall grid for the purposes of snowfall separation and bias adjustment in the frame of the present study.

Calibration of bias adjustment

The analysis of total precipitation ratios (RCM simulations with respect to gridded observations) for the EVAL period, which are computed to carry out the first step of the bias adjustment procedure, reveals substantial elevation dependencies. All simulations tend to overestimate total precipitation at high elevations (Fig. S3). This fact might ultimately be connected to an overestimation of surface snow amount in several EURO-CORDEX RCMs as reported by Terzago et al. (2017). As the precipitation ratio between simulations and observations depends approximately linearly on elevation, the calculation of $P_{\text{AF}}$ via a linear regression of the ratios against elevation (see Sect. 2.6) seems reasonable. By taking the inverse of this linear relation, $P_{\text{AF}}$ for every model and elevation can be derived. For the CCLM simulations, these correction factors do not vary much with height, while $P_{\text{AF}}$ for MPI-ESM–REMO and EC-EARTH–HIRHAM is much larger than 1 in low-lying areas, indicating a substantial underestimation of observed precipitation sums (Fig. 4a). However, for most elevations and simulations, $P_{\text{AF}}$ is generally smaller than 1, i.e. total precipitation is overestimated by the models. Similar model biases in the winter and spring seasons have already been reported in previous works (e.g. Rajczak et al., 2017; Smiatek et al., 2016). Especially at high elevations, these apparent positive precipitation biases could be related to observational undercatch, i.e. an underestimation of true precipitation sums by the observational analysis. Frei et al. (2003) estimated seasonal Alpine precipitation undercatch for three elevation intervals. Results show that measurement biases are largest in winter and increase with altitude. However, a potential undercatch (with a maximum of around 40 % at high elevations in winter; Frei et al., 2003) can only partly explain the often substantial overestimation of precipitation found in the present work.

Evaluation of snowfall indices in the EVAL period 1971–2005 for the 12 snowfall separated–bias-adjusted (${\text{RCM}}_{\text{sep}+\text{ba}}$) and 12 snowfall separated–not-bias-adjusted (${\text{RCM}}_{\text{sep}+\text{nba}}$) RCM simulations of the reduced model set vs. observation-based reference. Column (a) shows the mean September–May snowfall index statistics vs. elevation while the monthly snowfall indices (spatially averaged over the elevation intervals $<$ 1000, 1000–2000 and $>\mathrm{2000}$ $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) are displayed in columns (b) to (d).

[Figure omitted. See PDF]

After applying $P_{\text{AF}}$ to the daily precipitation fields, a snowfall fractionation at the initial $T^{\ast }$ of 2 ${}^{\circ }$C (see Eq. 4) would lead to a snowfall excess in all 12 simulations as models typically experience a cold winter temperature bias. To match the observation-based and spatio-temporally averaged reference snowfall below 2750 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, $T^{\ast }$ for all models needs to be decreased during the second step of the bias adjustment (Fig. 4b). The adjusted $T_{\mathrm{a}}^{\ast }$ values indicate a clear positive relation with the mean temperature bias in the EVAL period. This feature is expected since the stronger a particular model's cold bias is, the stronger the required adjustment of the snow fractionation temperature $T^{\ast }$ towards lower values is in order to avoid a positive snowfall bias. Various reasons for the scatter around a simple linear relation in Fig. 4b can be thought of. These include remaining spatial inaccuracies of the corrected precipitation grid, elevation-dependent temperature biases and misrepresented temperature–precipitation relationships on a daily scale. Note that precipitation and temperature biases heavily depend on the GCM–RCM chain and seem to be rather independent from each other. Concerning the partly substantial temperature biases of the EURO-CORDEX models shown in Fig. 4b, their magnitude largely agrees with Kotlarski et al. (2014; in reanalysis-driven simulations) and Smiatek et al. (2016).

Evaluation of snowfall indices

We next assess the performance of the bias adjustment procedure by comparing snowfall indices derived from separated–bias-adjusted RCM snowfall amounts against the observation-based reference. The period for which this comparison is carried out is EVAL; i.e. it is identical to the calibration period of the bias adjustment. We hence do not intend a classical cross validation exercise with separate calibration and validation periods, but we try to answer the following two questions. (a) Which aspects of the Alpine snowfall climate are adjusted? (b) For which aspects do biases remain even after application of the bias adjustment procedure?

Figure 5 shows the evaluation results of the six snowfall indices based on the separated–not-bias-adjusted simulated snowfall (${\text{RCM}}_{\text{sep}+\text{nba}}$) and the separated–bias-adjusted simulated snowfall (${\text{RCM}}_{\text{sep}+\text{ba}}$). In the first case the snowfall separation of raw precipitation is performed with $T^{\ast }=\mathrm{2}$ ${}^{\circ }$C, while in the second case precipitation is adjusted and the separation is performed with a bias-adjusted temperature $T_{\mathrm{a}}^{\ast }$. The first column (a, Fig. 5a) represents the mean September to May statistics, while columns (b) to (d) depict the seasonal cycle at monthly resolution for three distinct elevation intervals.

The analysis of $S_{\text{mean}}$ confirms that ${\text{RCM}}_{\text{sep}+\text{ba}}$ is able to reproduce the observation-based reference in the domain mean as well as in most individual elevation intervals. The domain-mean agreement is a direct consequence of the design of the bias adjustment procedure (see above). ${\text{RCM}}_{\text{sep}+\text{nba}}$, on the other hand, consistently overestimates $S_{\text{mean}}$ by up to a factor of 2.5 as a consequence of positive precipitation and negative temperature biases (see Fig. 4). Also, the seasonal cycle of $S_{\text{mean}}$ for ${\text{RCM}}_{\text{sep}+\text{ba}}$ yields a satisfying performance across all three elevation intervals, while ${\text{RCM}}_{\text{sep}+\text{nba}}$ tends to produce too much snowfall over all months and reveals an increasing model spread with elevation.

For the full domain and elevations around 1000 $\mathrm{m}$, the observation-based reference indicates a mean $S_{\text{freq}}$ of 20 % between September and May. Up to 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ ${\text{RCM}}_{\text{sep}+\text{ba}}$ reflects the increase in this index with elevation adequately. However, towards higher elevations the approximately constant $S_{\text{freq}}$ of 30 % in the reference is not captured by the simulation-derived snowfall. Notably during wintertime, both ${\text{RCM}}_{\text{sep}+\text{ba}}$ and ${\text{RCM}}_{\text{sep}+\text{nba}}$ produce too many snowfall days, i.e. overestimate snowfall frequency. This feature is related to the fact that climate models typically tend to overestimate the wet-day frequency over the Alps especially in wintertime (Rajczak et al., 2013) and that the bias adjustment procedure employed does not explicitly correct for potential biases in precipitation frequency. Due to the link between mean snowfall on one side and snowfall frequency and mean intensity on the other side, opposite results are obtained for the mean snowfall intensity $S_{\text{int}}$. ${\text{RCM}}_{\text{sep}+\text{ba}}$ largely underestimates mean intensities during snowfall days while ${\text{RCM}}_{\text{sep}+\text{nba}}$ typically better reflects the reference. Nevertheless, deviations during winter months at mid-elevations are not negligible. Mean September–May $S_{\text{frac}}$ in the reference exponentially increases with elevation. This behaviour is reproduced by both ${\text{RCM}}_{\text{sep}+\text{ba}}$ and ${\text{RCM}}_{\text{sep}+\text{nba}}$. Notwithstanding, ${\text{RCM}}_{\text{sep}+\text{ba}}$ results are more accurate compared to ${\text{RCM}}_{\text{sep}+\text{nba}}$, which turns out to be biased towards too-large snowfall fractions.

For the two heavy snowfall indices $S_{\mathrm{q}\mathrm{99}}$ and $S_{\mathrm{1}\mathrm{d}}$, ${\text{RCM}}_{\text{sep}+\text{nba}}$ appears to typically match the reference better than ${\text{RCM}}_{\text{sep}+\text{ba}}$. Especially at high elevations, ${\text{RCM}}_{\text{sep}+\text{ba}}$ produces too-low snowfall amounts. This again illustrates the fact that the bias adjustment procedure is designed to adjust biases in mean snowfall, but does not necessarily improve further aspects of the simulated snowfall climate.

The spatial patterns of $S_{\text{mean}}$ for the 12 RCM${}_{\text{sep}+\text{ba}}$ simulations from September to May are presented in Fig. 6. The observational-based reference (lower-right panel) reveals a snowfall distribution with highest values along the Alpine main ridge, whereas the Swiss plateau, southern Ticino and main valleys such as the Rhône and Rhine valley experience less snowfall. Almost all bias-adjusted models are able to represent the overall picture with snow-poor lowlands and snow-rich Alpine regions. Nevertheless, substantial differences in the observations concerning the spatial snowfall pattern can arise. EC-EARTH–HIRHAM, for example, is subject to a noisy structure. This could be the result of frequent grid-cell storms connected to parameterisations struggling with complex topographies. Such inaccuracies in the spatial pattern are not corrected for by our simple bias adjustment approach which only targets domain-mean snowfall amounts at elevations below 2750 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ and that does not considerably modify the simulated spatial snowfall patterns. Note that these patterns are obviously strongly determined by the RCM itself and only slightly depend on the driving GCM (see, for instance, the good agreement among the CCLM and the RCA simulations).

Spatial distribution of mean September–May snowfall, $S_{\text{mean}}$, in the EVAL period 1971–2005 and for the 14 snowfall separated–bias-adjusted RCM simulations (${\text{RCM}}_{\text{sep}+\text{ba}}$) of the full model set. Lower-right panel: observation-based reference.

[Figure omitted. See PDF]

In summary, after applying the bias adjustment to the simulations, most snowfall indices are fairly well represented at elevations below 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ With increasing altitude and smaller sample sizes in terms of number of grid cells, reference and ${\text{RCM}}_{\text{sep}+\text{ba}}$ diverge. This might be caused by the remaining simulated overestimation of $S_{\text{freq}}$ and an underestimation of $S_{\text{int}}$. While the bias adjustment approach leads to a reduction of $S_{\text{int}}$ due to the total precipitation adjustment, $S_{\text{freq}}$ is only slightly modified by this correction and by the adjustment of $T^{\ast }$. Nevertheless, these two parameters strongly influence other snowfall indices. The counteracting effects of overestimated $S_{\text{freq}}$ and underestimated $S_{\text{int}}$ result in appropriate amounts of $S_{\text{mean}}$, whereas discrepancies for $S_{\mathrm{q}\mathrm{99}}$ and $S_{\mathrm{1}\mathrm{d}}$ are mainly driven by the underestimation of $S_{\text{int}}$.

Snowfall projections for the late 21st century

For the study of climate change signals, the analysis domain is extended to the entire Alps (see Sect. 2.3). Due to the identified difficulties of bias-adjusting certain snowfall indices (see Sect. 3.4), the emphasis is on the relative signals of change (see Eq. 2). This type of change can be expected to be less dependent on the remaining inaccuracies after the adjustment. If not stated otherwise, all results in this section are based on the ${\text{RCM}}_{\text{sep}+\text{ba}}$ data, i.e. on separated–bias-adjusted RCM snowfall, and on the RCP8.5 emission scenario. Depending on the type of analysis, either the full or the reduced model set is used (see Table 1 and Supplement, Part B).

Projections for seasonal $S_{\text{mean}}$ show a considerable decrease over the entire Alpine domain (Fig. 7). Most RCMs project largest percentage losses of more than 80 % across the Alpine forelands and especially in its topographic depressions such as the Po and Rhône valleys. Over the Alpine ridge, reductions are smaller but still mostly negative. Elevated regions between southeastern Switzerland, northern Italy and Austria seem to be least affected by the overall snowfall reduction. Some of the simulations (e.g. CNRM–RCA, MPI-ESM–RCA or MPI-ESM–REMO) project only minor changes in these regions. Experiments employing the same RCM but different driving GCMs (e.g. the four simulations of RCA) but also experiments employing the same GCM but different RCMs (e.g. the four simulations driven by EC-EARTH, though different realisations) can significantly disagree in regional-scale change patterns and especially in the general magnitude of change. This highlights a strong influence of both the driving GCMs and the RCMs themselves on snowfall changes, representing effects of large-scale circulation and mesoscale response, respectively.

Spatial distribution of relative changes (SCEN period 2070–2099 with respect to CTRL period 1981–2010) in mean September–May snowfall, $\mathit{\delta }{S}_{\text{mean}}$, for RCP8.5 and for the 14 snowfall separated–bias-adjusted RCM simulations (${\text{RCM}}_{\text{sep}+\text{ba}}$) of the full model set. For RCP4.5, see Fig. S7.

[Figure omitted. See PDF]

A more detailed analysis is provided in Fig. 8 which addresses the vertical and seasonal distribution of snowfall changes. It reveals that relative (seasonal mean) changes in $S_{\text{mean}}$ appear to be strongly dependent on elevation (Fig. 8, top left panel). The multi-model mean change ranges from $-$80 % at low elevations to $-$10 % above 3000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ The largest differences between neighbouring elevation intervals are obtained from 750 to 1500 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ Over the entire Alps, the results show a reduction of $S_{\text{mean}}$ by $-$35 to $-$55 % with a multi-model mean of $-$45 %. The multi-model spread appears to be rather independent of elevation and is comparably small, confirming that, overall, the spatial distribution of the change pattern is similar across all model chains (see Fig. 7). All simulations point to decreases over the entire 9-month period from September to May for the two elevation intervals $<$ 1000 and 1000 to 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ Above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, individual simulations show an increase in $S_{\text{mean}}$ by up to 20 % in midwinter which leads to a slightly positive change in the multi-model mean in January and February.

Relative changes (SCEN period 2070–2099 with respect to CTRL period 1981–2010) in snowfall indices based on the 14 snowfall separated–bias-adjusted RCM simulations (${\text{RCM}}_{\text{sep}+\text{ba}}$) of the full model set for RCP8.5. Column (a) shows the mean September–May snowfall index statistics vs. elevation. Monthly snowfall index changes (spatially averaged over the elevation intervals $<\mathrm{1000}$, 1000–2000 and $>\mathrm{2000}$ $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) are displayed in columns (b) to (d).

[Figure omitted. See PDF]

Decreases in $S_{\text{freq}}$ are very similar to changes in mean snowfall. Mean September–May changes are largest below 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, while differences among elevation intervals become smaller at higher elevations. In-between is a transition zone with rather strong changes with elevation, which approximately corresponds to the mean elevation of the September–May zero-degree line in today's climate (e.g. Ceppi et al., 2012; MeteoSwiss, 2016). Individual simulations with large reductions in $S_{\text{mean}}$, such as the RCA experiments, also project the strongest declines in $S_{\text{freq}}$. In contrast, the mean snowfall intensity $S_{\text{int}}$ is subject to the smallest percentage variations in our set of snowfall indices. Strong percentage changes for some models in September are due to the small sample size (only few grid points considered) and the low snowfall amounts in this month. Apart from mid-elevations with decreases of roughly $-$10 %, mean intensities from September to May are projected to remain almost unchanged by the end of the century. For both seasonal and monthly changes, model agreement is best for high elevations while the multi-model spread is largest for the lowlands. The large model spread at low elevations might be caused by the small number of grid points used for averaging over the respective elevation interval, especially in autumn and spring.

Intercomparison of various snowfall indices and relationship with monthly mean temperature in CTRL. For each panel, the monthly mean statistics for each 250 $\mathrm{m}$ elevation interval and for each of the 12 snowfall separated–bias-adjusted (${\text{RCM}}_{\text{sep}+\text{ba}}$) RCM simulations of the reduced model set (${\text{RCM}}_{\text{sep}+\text{ba}}$) were derived (black circles). Red triangles denote the multi-model mean for a specific month and elevation interval. The monthly statistics were calculated by considering all grid points of the specific elevation intervals which are available for both variables in the corresponding scatter plot only (area consistency). Relative changes are based on the RCP8.5-driven simulations (SCEN 2070–2099 with regard to CTRL 1981–2010).

[Figure omitted. See PDF]

Comparison of temperature probability, number of snowfall days, mean snowfall intensity and snowfall frequency for the CTRL period 1981–2010 and SCEN period 2070–2099 for RCP8.5. The analysis is based on data from the 12 snowfall separated–bias-adjusted RCM simulations (${\text{RCM}}_{\text{sep}+\text{ba}}$) of the reduced model set. The top row (a) depicts the PDF of the daily temperature distribution, while the second row (b) shows the mean number of snowfall days between September and May, i.e. days with $S>\mathrm{1}$ $\mathrm{mm}$ (see Table 2), in a particular temperature interval. The third row (c) represents the mean snowfall intensity, $S_{\text{int}}$, for a given snowfall temperature interval. In addition, the Clausius–Clapeyron relationship, centred at the $-$10 ${}^{\circ }$C mean $S_{\text{int}}$ for SCEN, is displayed by the black dashed line. Row (d) shows the snowfall frequency, $S_{\mathrm{freq}}$, for a given temperature interval. PDFs and mean $S_{\text{int}}$ and $S_{\mathrm{freq}}$ were calculated by creating daily mean temperature bins of width 1 ${}^{\circ }$C.

[Figure omitted. See PDF]

Similar results are obtained for the heavy snowfall indices $S_{\mathrm{q}\mathrm{99}}$ and $S_{\mathrm{1}\mathrm{d}}$. While percentage decreases at lowermost elevations are even larger than for $S_{\text{mean}}$, losses at high elevations are less pronounced, resulting in similar domain-mean change signals for heavy and mean snowfall. Substantial differences between monthly $\mathit{\delta }{S}_{\mathrm{q}\mathrm{99}}$ and $\mathit{\delta }{S}_{\mathrm{1}\mathrm{d}}$ appear at elevations below 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ Here, percentage losses of $S_{\mathrm{q}\mathrm{99}}$ are typically slightly more pronounced. Above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ both indices appear to remain almost constant between January and March with change signals close to zero. The multi-model mean changes even hint to slight increases in both indices. Concerning changes in the snowfall fraction, i.e. in the relative contribution of snowfall to total precipitation, our results indicate that current seasonal and domain-mean $S_{\text{frac}}$ might drop by about $-$50 % (Fig. 8, lowermost row). Below 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, the strength of the signal is almost independent of the month, and multi-model average changes in the snow fraction of about $-$80 % are obtained. At higher elevations changes during midwinter are less pronounced compared to autumn and spring but still negative.

Discussion

Effect of temperature, snowfall frequency and intensity on snowfall changes

The results in Sect. 4 indicate substantial changes in snowfall indices over the Alps in regional climate projections. With complementary analyses presented in Figs. 9 and 10, we shed more light on the responsible mechanisms, especially concerning projected changes in mean and heavy snowfall. For this purpose Fig. 9a, b and e–f show the relationship of both mean and heavy snowfall amounts in the CTRL period and their respective percentage changes with the climatological CTRL temperature of the respective (climatological) month, elevation interval and GCM–RCM chain. For absolute amounts ($S_{\text{mean}}$, $S_{\mathrm{q}\mathrm{99}}$; Fig. 9a and e) a clear negative relation is found, i.e. the higher the CTRL temperature, the lower the snowfall amounts. For $S_{\text{mean}}$ the relation levels off at mean temperatures higher than about 6 ${}^{\circ }$C with mean snowfall amounts close to zero. For temperatures below about $-$6 ${}^{\circ }$C a considerable spread of snowfall amounts is obtained, i.e. mean temperature does not seem to be the controlling factor here. Relative changes in both quantities (Fig. 9b and f), however, are strongly controlled by the CTRL period's temperature level with losses close to 100 % for warm climatic settings and partly increasing snowfall amounts for colder climates. This dependency of relative snowfall changes on CTRL temperature is in line with previous works addressing future snowfall changes on both hemispheric and regional scales (de Vries et al., 2014; Krasting et al., 2013; Räisänen, 2016). The spread of changes within a given CTRL temperature bin can presumably be explained by the respective warming magnitudes that differ between elevations, months and GCM–RCM chains. About half of this spread can be attributed to the month and the elevation alone (compare the spread of the black markers to the one of the red markers which indicate multi-model averages).

For most months and elevation intervals, percentage reductions in $S_{\text{mean}}$ and $S_{\mathrm{q}\mathrm{99}}$ reveal an almost linear relationship with $\mathit{\delta }{S}_{\text{freq}}$ (Fig. 9c and g). The decrease in $S_{\text{freq}}$ with future warming can be explained by a shift in the temperature probability distribution towards higher temperatures, leading to fewer days below the freezing level (Fig. 10a). Across the three elevation intervals $<\mathrm{1000}$, 1000–2000 and $>\mathrm{2000}$ $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, relative changes in the number of days with temperatures below the freezing level ($T\le \mathrm{0}$ ${}^{\circ }$C) are in the order of $-$65, $-$40 and $-$20 %, respectively (not shown). This approximately corresponds to the simulated decrease in $S_{\text{freq}}$ (see Fig. 8), which in turn is of a similar magnitude as that found in previous works addressing future snowfall changes in the Alps (Schmucki et al., 2015b; Zubler et al., 2014). Due to the general shift in the temperature distribution and the loss of very cold days (Fig. 10a), future snowfall furthermore occurs in a narrower temperature range (Fig. 10b).

Contrasting this general pattern of frequency-driven decreases in both mean and heavy snowfall, no changes or even slight increases in $S_{\text{mean}}$, $S_{\mathrm{q}\mathrm{99}}$ and $S_{\mathrm{1}\mathrm{d}}$ at high elevations are expected in midwinter (see Fig. 8). This can to some part be explained by the general increase in total winter precipitation (Rajczak et al., 2017; Smiatek et al., 2016) that obviously offsets the warming effect in high-elevation regions where a substantial fraction of the future temperature PDF is still located below the rain–snow transition (Fig. 10, top row). This process has also been identified in previous works to be, at last partly, responsible for future snowfall increases (de Vries et al., 2014; Krasting et al., 2013; Räisänen, 2016). Furthermore, the magnitude of the increases in both mean and heavy snowfall is obviously driven by positive changes in $S_{\text{int}}$, while $S_{\text{freq}}$ remains constant (Fig. 9c and g). An almost linear relationship between positive changes in $S_{\text{int}}$ and positive changes in $S_{\text{mean}}$ and $S_{\mathrm{q}\mathrm{99}}$ is obtained (Fig. 9d and h; upper-right quadrants. Nevertheless, the high-elevation midwinter growth in $S_{\text{mean}}$ is smaller than the identified increases in mean winter total precipitation. This can be explained by the persistent decrease in $S_{\text{frac}}$ during the cold season (see Fig. 8, lowermost row).

For elevation intervals with simulated monthly temperatures between $-$6 and 0 ${}^{\circ }$C in the CTRL period, $S_{\text{mean}}$ appears to decrease stronger than $S_{\mathrm{q}\mathrm{99}}$ (see Fig. 9b and f). O'Gorman (2014) found a very similar behaviour when analysing mean and extreme snowfall projections over the Northern Hemisphere within a set of GCMs. This finding is related to the fact that future snowfall decreases are mainly governed by a decrease in snowfall frequency while snowfall increases in high-elevated regions in midwinter seem to be caused by increases in snowfall intensity. This can obviously be explained by the insensitivity of the temperature interval at which extreme snowfall occurs to climate warming and by the shape of the temperature–snowfall-intensity distribution itself (Fig. 10c). The likely reason behind the positive changes in $S_{\text{int}}$ at high-elevated and cold regions is the higher water holding capacity of the atmosphere in a warmer climate. According to the Clausius–Clapeyron relation, saturation vapour pressure increases by about 7 $\mathrm{\%}{\mathrm{degree}}^{-\mathrm{1}}$ warming (Held and Soden, 2006). Previous studies have shown that simulated changes in heavy and extreme precipitation on daily timescales are consistent with this theory (e.g. Allen and Ingram, 2002; Rajczak et al., 2017). In terms of snowfall, we find the Clausius–Clapeyron relation to be applicable for negative temperatures up to approximately $-$5 ${}^{\circ }$C as well (Fig. 10c, dashed lines). Inconsistencies for temperatures between $-$5 and 0 ${}^{\circ }$C are due to a snow fraction $f_{\mathrm{s}}$ $<$ 100 % for corresponding precipitation events.

Schematic illustration of the control of changes in snowfall intensity on changes in mean and extreme snowfall. (a) Relation between temperature at snowfall and mean snowfall intensity. (b) Daily temperature PDF for a warm control climate (low elevations or transition seasons, i.e. beginning or end of winter). (c) Daily temperature PDF for a cold control climate (high elevations or midwinter). The blue line denotes the historical CTRL period; the red line denotes the future SCEN period. The light grey shaded area represents the overall temperature interval in which snowfall occurs; the dark grey shading shows the preferred temperature interval for heavy snowfall to occur.

[Figure omitted. See PDF]

Similar to Fig. 8 but showing projected changes in mean snowfall, $\mathit{\delta }{S}_{\text{mean}}$, and heavy snowfall, $\mathit{\delta }{S}_{\mathrm{q}\mathrm{99}}$, for the emission scenarios RCP4.5 and 8.5 and based on the reduced model set. See Fig. S10 for the emission scenario uncertainty of the remaining four snowfall indices. See Fig. S15 for the respective figure based on the full model set.

[Figure omitted. See PDF]

For further clarification, Fig. 11 schematically illustrates the governing processes behind the changes in mean and heavy snowfall that differ between climatologically warm (decreasing snowfall) and climatologically cold climates (increasing snowfall). As shown in Fig. 10 (third row), the mean $S_{\text{int}}$ distribution is rather independent of future warming and similar temperatures are associated with similar mean snowfall intensities. In particular, the heaviest snowfall is expected to occur slightly below the freezing level in both the CTRL and the SCEN period (Fig. 11a). How often do such conditions prevail in the two periods? In a warm current climate, i.e. at low elevations or in the transition seasons, heavy snowfall only rarely occurs as the temperature interval for highest snowfall intensity is already situated in the left tail of the CTRL period's temperature distribution (Fig. 11b). With future warming, i.e. with a shift in the temperature distribution to the right, the probability for days to occur in the heavy snowfall temperature interval (dark grey shading) decreases more strongly than the probability for days to occur in the overall snowfall regime (light grey shading). This results in (1) a general decrease in snowfall frequency, (2) a general decrease in mean snowfall intensity and (3) a general and similar decrease in both mean and heavy snowfall amounts. In contrast, at cold and high-elevated sites, CTRL period temperatures are often too low to trigger heavy snowfall since a substantial fraction of the temperature PDF is located to the left of the heavy snowfall temperature interval (Fig. 11c). The shifted distribution in a warmer SCEN climate, however, peaks within the temperature interval that favours heavy snowfall. This leads to a probability increase for days to occur in the heavy snowfall temperature range despite the general reduction in $S_{\text{freq}}$ (lower overall probability for days to occur in the entire snowfall regime, light grey). As a consequence, mean $S_{\text{int}}$ tends to increase and the reduction of heavy snowfall amounts is less pronounced (or even of opposing sign) than the reduction in mean snowfall. For individual (climatologically cold) regions and seasons, the increase in mean $S_{\text{int}}$ might even compensate the $S_{\text{freq}}$ decrease, resulting in an increase in both mean and heavy snowfall amounts. Note that in a strict sense these explanations only hold in the case that the probability for snowfall to occur at a given temperature does not change considerably between the CTRL and the SCEN period. This behaviour is approximately found (Fig. 10d), which presumably indicates only minor contributions of large-scale circulation changes and associated humidity changes on both the temperature–snowfall frequency and the temperature–snowfall intensity relation.

Emission scenario uncertainty

The projections presented in the previous sections are based on the RCP8.5 emission scenario, but will depend on the specific scenario considered. To assess this uncertainty we here compare the ${\text{RCM}}_{\text{sep}+\text{ba}}$ simulations for the previously shown RCP8.5 emission scenario against those assuming the more moderate RCP4.5 scenario. As a general picture, the weaker RCP4.5 scenario is associated with less pronounced changes in snowfall indices (Fig. 12). Differences in mean seasonal $\mathit{\delta }{S}_{\text{mean}}$ between the two emission scenarios are most pronounced below 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ where percentage changes for RCP4.5 are about one third smaller than for RCP8.5. At higher elevations, multi-model mean changes agree better and the multi-model ranges for the two emission scenarios start overlapping; i.e. individual RCP4.5 experiments can be located in the RCP8.5 multi-model range and vice versa. Over the entire Alpine domain, about $-$25 % of current snowfall is expected to be lost under the moderate RCP4.5 emission scenario while a reduction of approximately $-$45 % is projected for RCP8.5. For seasonal cycles, the difference of $\mathit{\delta }{S}_{\text{mean}}$ between RCP4.5 and RCP8.5 is similar for most months and slightly decreases with altitude. Above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, the simulated increase in $S_{\text{mean}}$ appears to be independent of the chosen RCP in January and February, while negative changes before and after midwinter are more pronounced for RCP8.5. Alpine domain-mean $\mathit{\delta }{S}_{\mathrm{q}\mathrm{99}}$ almost doubles under the assumption of stronger GHG emissions. This is mainly due to differences at low elevations, whereas above 2000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ $\mathit{\delta }{S}_{\mathrm{q}\mathrm{99}}$ does not seem to be strongly affected by the choice of the emission scenario. Differences in monthly mean changes are in close analogy to $\mathit{\delta }{S}_{\text{mean}}$. Higher emissions lead to a further negative shift in $\mathit{\delta }{S}_{\mathrm{q}\mathrm{99}}$. Up to mid-elevations differences are rather independent of the season. However, at highest elevations and from January to March, differences between RCP4.5 and RCP8.5 are very small.

Despite the close agreement of midwinter snowfall increases at high elevations between the two emission scenarios, obvious differences in the spatial extent of the region of mean seasonal snowfall increases can be found (see Figs. S7 and 7 for $\mathit{\delta }{S}_{\text{mean}}$, and Figs. S8 and S9 for $\mathit{\delta }{S}_{\mathrm{q}\mathrm{99}}$). In most simulations, the number of grid cells along the main Alpine ridge that show either little change or even increases in seasonal mean $S_{\text{mean}}$ or $S_{\mathrm{q}\mathrm{99}}$ is larger for RCP4.5 than for RCP8.5 with its larger warming magnitude.

Relative and absolute changes (SCEN period 2070–2099 with respect to CTRL period 1981–2010) in mean September–May snowfall indices based on a subset of seven snowfall separated–bias-adjusted (${\text{RCM}}_{\text{sep}+\text{ba}}$), seven snowfall separated–not-bias-adjusted (${\text{RCM}}_{\text{sep}+\text{nba}}$) and seven raw snowfall RCM simulations (${\text{RCM}}_{\text{raw}}$) for RCP8.5. Only RCM simulations of the reduced model set providing raw snowfall as an output variable (see Table 1) were used in this analysis.

[Figure omitted. See PDF]

Intercomparison of projections with separated and raw snowfall

The snowfall projections presented above are based on the ${\text{RCM}}_{\text{sep}+\text{ba}}$ data set, i.e. on separated–bias-adjusted snowfall amounts. To assess the robustness of these estimates we here compare the obtained change signals against the respective signals based on ${\text{RCM}}_{\text{sep}+\text{nba}}$ (separated–not-bias-adjusted) and simulated raw snowfall output (${\text{RCM}}_{\text{raw}}$). This comparison is restricted to the seven RCMs providing raw snowfall as an output variable (see Table 1).

The three different change estimates agree well with each other in terms of relative snowfall change signals (Fig. 13, top row). Multi-model mean relative changes are very similar for all analysed snowfall indices and elevation intervals. In many cases, ${\text{RCM}}_{\text{sep}+\text{nba}}$ is subject to slightly smaller percentage decreases. Multi-model mean differences between ${\text{RCM}}_{\text{sep}+\text{ba}}$, ${\text{RCM}}_{\text{sep}+\text{nba}}$ and ${\text{RCM}}_{\text{raw}}$ simulations are smaller than the corresponding multi-model spread of ${\text{RCM}}_{\text{sep}+\text{ba}}$ simulations and emission scenario uncertainties (see Figs. 12, 13 and S10). This agreement in terms of relative change signals is in contrast to absolute change characteristics (Fig. 13, bottom row). Results based on the three data sets agree in the sign of change, but not in their magnitude, especially at high elevations $>\mathrm{2000}$ $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ As the relative changes are almost identical, the absolute changes strongly depend upon the treatment of biases in the control climate.

In summary, these findings indicate that (a) the snowfall separation method developed in the present work yields rather good proxies for relative changes in snowfall indices in raw RCM output (which is not available for all GCM–RCM chains) and that (b) the additional bias adjustment of separated snowfall amounts only has a weak influence on relative change signals of snowfall indices but can have substantial effects on absolute changes.

Conclusions and outlook

The present work makes use of state-of-the-art EURO-CORDEX RCM simulations to assess changes in snowfall indices over the European Alps by the end of the 21st century. For this purpose, snowfall is separated from total precipitation using near-surface air temperature in both the RCMs and in the observation-based estimate on a daily basis. The analysis yields a number of robust signals, consistent across a range of climate model chains and across emission scenarios. Relating to the main objectives we find the following:

Snowfall separation on the RCM grid. Binary snow fractionation with a fixed temperature threshold on coarse-resolution grids (with 12 $\mathrm{km}$ resolution) leads to an underestimation of mean snowfall and an overestimation of heavy snowfall. To overcome these deficiencies, the Richards snow fractionation method is implemented. This approach expresses that the coarse-grid snow fraction depends not only on daily mean temperature, but also on topographical subgrid-scale variations. Accounting for the latter results in better estimates for mean and heavy snowfall. However, due to limited observational coverage the parameters of this method are fitted for Switzerland only and are then applied to the entire Alpine domain. Whether this spatial transfer is robust could be further investigated by using observational data sets covering the full domain of interest but is out of the scope of this study.

Snowfall bias adjustment. Simulations of the current EURO-CORDEX ensemble are subject to considerable biases in precipitation and temperature, which translate into biased snowfall amounts. In the EVAL period, simulated precipitation is largely overestimated, with increasing biases toward higher altitudes. On the other hand, simulated near-surface temperatures are generally too low with the largest deviations over mountainous regions. These findings were already reported in previous studies for both the current EURO-CORDEX data set and for previous RCM ensembles (e.g. Frei et al., 2003; Kotlarski et al., 2012, 2015; Rajczak et al., 2013; Smiatek et al., 2016). By implementing a simple bias adjustment approach, we are able to partly reduce these biases and the associated model spread, which should enable more robust change estimates. The adjusted model results reproduce the seasonal cycles of mean snowfall fairly well. However, substantial biases remain in terms of heavy snowfall, snowfall intensities (which in general are overestimated), snowfall frequencies and spatial snowfall distributions. Further improvements might be feasible by using more sophisticated bias adjustment methods, such as quantile mapping (e.g. Rajczak et al., 2016), local intensity scaling of precipitation (e.g. Schmidli et al., 2006) or weather generators (e.g. Keller at al., 2016). The advantages of the approach employed here are its simplicity, its direct linkage to the snowfall separation method and, as a consequence, its potential ability to account for non-stationary snowfall biases. Furthermore, a comparison to simulated raw snowfall for a subset of seven simulations revealed that relative change signals are almost independent of the chosen post-processing strategy.

Snowfall projections for the late 21st century. Snowfall climate change signals are assessed by deriving the changes in snowfall indices between the CTRL period 1981–2010 and the SCEN period 2070–2099. Our results show that by the end of the 21st century, snowfall over the Alps will be considerably reduced. Between September and May mean snowfall is expected to decrease by approximately $-$45 % (multi-model mean) under an RCP8.5 emission scenario. For the more moderate RCP4.5 scenario, multi-model mean projections show a decline of $-$25 %. These results are in good agreement with previous works (e.g. de Vries et al., 2014; Piazza et al., 2014; Räisänen, 2016). Low-lying areas experience the largest percentage changes of more than $-$80 %, while the highest Alpine regions are only weakly affected. Variations of heavy snowfall, defined by the 99 % all-day-snowfall percentile, show an even more pronounced signal at low-lying elevations. With increasing elevation, percentage changes in heavy snowfall are generally smaller than for mean snowfall. O'Gorman (2014) found a very similar behaviour by analysing projected changes in mean and extreme snowfall over the entire Northern Hemisphere. He pointed out that heavy and extreme snowfall occurs near an optimal temperature (near or below freezing, but not too cold), which seems to be independent of climate warming. We here confirm this finding. At mid- and high elevations, heavy snowfall in a warmer climate will still occur in the optimal temperature range; hence, heavy snowfall amounts will decrease less strongly compared to mean snowfall and may even increase in some areas.

At first approximation, the magnitude of future warming strongly influences the reduction of mean and heavy snowfall by modifying the snowfall frequency. Snowfall increases may occur at high (and thus cold) elevations, and these are not caused by frequency changes. Here, snowfall increases due to (a) a general increase in total winter precipitation combined with only minor changes in snowfall frequency and (b) more intense snowfall. This effect has a pronounced altitudinal distribution and may be particularly strong under conditions (depending upon location and season) where the current climate is well below freezing. With the expected warming, a shift towards a temperature range more favourable to snowfall (near or below freezing, but not too cold) can be expected with corresponding increases in mean snowfall, despite a general decrease in the snowfall fraction.

Note that individual EURO-CORDEX experiments were completely or partly omitted from our analysis (see Supplement, Part B, and Table 1). However, when comparing the ensemble results based on the reduced and the full model set, only slight differences are found and our main results and conclusions do not change (compare Fig. 12 to Fig. S15). This indicates a robust ensemble analysis that is only slightly affected by potentially critical shortcomings of individual simulations.

The identified future changes in snowfall over the Alps can lead to a variety of impacts in different sectors. With decreasing snowfall frequencies and the general increase in the snowline (e.g. Beniston, 2003; Gobiet et al., 2014; Hantel et al., 2012), both associated with temperature changes, ski lift operators are looking at an uncertain future. A shorter snowfall season will likely put them under greater financial pressure. Climate change effects might be manageable only for ski areas reaching up to high elevations (e.g. Elsasser and Bürki, 2002). Even so, these resorts might start operating later in the ski season, and the snow conditions into early spring could change less dramatically due to projected high-elevation snowfall increases in midwinter. A positive aspect of the projected decrease in snowfall frequency might be reduced expenditures for airport and road safety (e.g. Zubler et al., 2015).

At lower altitudes, an intensification of winter precipitation, combined with smaller snowfall fractions (Serquet et al., 2013), increases the flood potential (Beniston, 2012). Snow can act as a buffer by releasing melt water constantly over a longer period of time. With climate warming, this storage capacity is lost, and heavy precipitation immediately drains into streams and rivers which might not be able to take up the vast amount of water fast enough. Less snowmelt will also have impacts on hydropower generation and water management (e.g. Weingartner et al., 2013). So far, many Alpine regions have been able to bypass dry periods by tapping into melt water from mountainous regions. With reduced snowpacks due to less snowfall, water shortage might become a serious problem in some areas.

Regarding specific socio-economic impacts caused by extreme snowfall events, conclusions based on the results presented in this study are difficult to draw. It might be possible that the 99 % all-day-snowfall percentile we used for defining heavy snowfalls is not appropriate for speculating about future evolutions of (very) rare events (Schär et al., 2016). To do so, one might consider applying a generalised extreme value analysis which is more suitable for answering questions related to rare extreme events.

The EURO-CORDEX RCM data analysed in the present work are publicly available – parts of them for non-commercial use only – via the Earth System Grid Federation archive (ESGF; e.g. https://esgf-data.dkrz.de). The observational data sets RhiresD and TabsD as well as the snow depth data for Switzerland are available for research and educational purposes from [email protected]. The analysis code is available from the corresponding author on request.

The Supplement related to this article is available online at https://doi.org/10.5194/tc-12-1-2018-supplement.

The authors declare that they have no conflict of interest.

Acknowledgements

We gratefully acknowledge the support of Jan Rajczak, Urs Beyerle and Curdin Spirig (ETH Zurich) as well as Elias Zubler (MeteoSwiss) in data acquisition and preprocessing. Christoph Frei (MeteoSwiss) and Christoph Marty (WSL-SLF) provided important input on specific aspects of the analysis. The GTOPO30 digital elevation model is available from the US Geological Survey. Finally, we thank the climate modelling groups of the EURO-CORDEX initiative for producing and making their model output available. Edited by: Xavier Fettweis Reviewed by: three anonymous referees