On the effect of short-term climate variability on mountain glaciers: insights from a case study

Abstract

Studies addressing the response of glaciers to climate change have so far analyzed the effect of long-term trends in a particular set of meteorological variables only, implicitly assuming an unaltered climatic variability. Here a framework for distinguishing between year-to-year, month-to-month and day-to-day variability is proposed. Synthetically generated temperature and precipitation time series following the same long-term trend but with altered variability are then used to force an ice-dynamics model set up for Rhonegletscher, Swiss Alps. In the case of temperature, variations in the day-to-day variability are shown to have a larger effect than changes at the yearly scale, while in the case of precipitation, variability changes are assessed as having negligible impact. A first set of scenarios is used to show that compared to reference, doubling the temperature variability can reduce glacier ice volume by up to 64% within half a decade. A second set derived from the results of the European ENSEMBLES project, however, shows that such changes are expected to remain below 8% even for extreme scenarios. Although the latter results relativize the importance of the effect in the near future, the analyses indicate that at least caution is required when assuming ‘unchanged variability’.

Introduction

Glaciers are known to be robust natural climate indicators (e.g. Dyurgerov and Meier, 2000; Barry, 2006) and have often been used to infer past climate evolution at decadal to millennial timescales (e.g. Oerlemans, 1994, 2005). At the same time, a number of studies have assessed the impact of changing climate on glacier evolution on a global (e.g. Meier, 1984; Gregory and Oerlemans, 1998; Raper and Braithwaite, 2006), regional (e.g. Ohmura and others, 1996; Schneeberger and others, 2003; Huss, 2011) or local scale (e.g. Horton and others, 2006; Stahl and others, 2008; Farinotti and others, 2012), and the topic is still a research focus (e.g. Gabbi and others, 2012; Huss, 2012; Springer and others, 2013). Commonly, however, future glacier evolution is assessed only with respect to a change in the mean of a particular variable – temperature or precipitation in most cases – whereas changes in the variability are either accounted for implicitly by forcing a given glacier model with the output of some climate model directly (e.g. Matulla and others, 2009) or are simply neglected (e.g. Oerlemans and others, 1998; Wright and others, 2005; Aðalgeirsdóttir and others, 2011).

To date, only a few studies have addressed the role of climate variability in glacier evolution, and these have focused exclusively on fluctuations with periodicities larger than one season. Nye (1960) used kinematic-wave theory to analyze the response of a glacier to seasonal and climatic periodic changes in accumulation, and showed that the interference of the waves generating from different glacier parts gives wide scope for variation in glacier evolution. Focusing on the glaciers around Mount Baker, Cascades Mountains, USA, and using glacier modeling, Roe and O’Neal (2009) pointed out how natural variability alone is capable of producing kilometer-scale (~40%) excursions in glacier length on multi-decadal and centennial timescales. A similar conclusion was drawn earlier by Reichert and others (2002) who analyzed two glaciers in Norway and the Swiss Alps, respectively, suggesting that preindustrial fluctuations of the glaciers, including their advance during the Little Ice Age, can be explained by internal variability in the climate system alone. Both studies emphasized the difficulty of separating a climate-change signal from natural variability, a topic that was addressed in more detail by Pederson and others (2004) and Roe (2011), and that implicitly highlights the importance of correctly addressing the evolution of climate variability when aiming at future projections. More generally, the importance of properly addressing variability changes in climate-change impact studies, especially studies of extreme events, was pointed out by Katz and Brown (1992) who used statistical theory to demonstrate that the frequency of extreme events is relatively more dependent on changes in the variability than in the mean of the climate.

This paper examines the effect that various assumptions about short-term climate variability (i.e. variability at time-scales shorter than the glacier response time) have when carrying out projections for future glacier evolution. The effects of an altered year-to-year, month-to-month and day-to-day variability are analyzed separately for temperature and precipitation by forcing a three-dimensional, full-Stokes ice-dynamics model with ensembles of meteorological time series obtained with a weather generator developed ad hoc. Two sets of scenarios are derived. In the first, currently observed variability is doubled or halved, making it possible to demonstrate the large potential effects. The second set, derived from the results of the European ENSEMBLES project (Van der Linden and Mitchell, 2009), is then used to assess the likely effect within the next century. The analyses are carried out for Rhonegletscher, Swiss Alps, but may be generalized for other locations.

Study Site and Data

Rhonegletscher is a south-southwest-exposed, medium-sized valley glacier (area ~16 km² as of 2007) in the main ridge of the Swiss Alps (Fig. 1a and b). Its elevation range in 2007 was 3208–3597 m a.s.l., with an average equilibrium-line altitude for the period 1865–2006 of 2950 m a.s.l. (Huss and others, 2008). Due to its good accessibility, the glacier has been the object of numerous investigations, starting with ice-flow speed measurements by Mercanton (1916), and ranging from analyses of the glacier mass balance (e.g. Müller and others, 1980; Chen and Funk, 1990; Huss and others, 2008) or the hydrology of the basin (e.g. Kasser, 1973; Bernath, 1991; Klok and others, 2001; Farinotti and others, 2012) to studies of the ice dynamics (e.g. Haeberli and Schlüchter, 1987; Sugiyama and others, 2007; Jouvet and others, 2009, 2011a) and internal ice deformation (e.g. Sugiyama and others, 2008; Keller and Blatter, 2012).

Fig. 1. (a) Glacier outline and surface contours for Rhonegletscher as of 2007. The blue line indicates the longitudinal section shown in Figures 5c and d and 6c and d. The black cross is a spatial reference (World Geodetic System 1984 coordinates). (b) Overview of Rhonegletscher, and location inside Switzerland. The image refers to September 2007. (c) Mean annual evolution of daily air temperature and monthly precipitation for the reference period 1980–2009. The data refer to the glacier center point (46.61° N, 8.39° E; 2700 m a.s.l.).

The surface topography of Rhonegletscher is known for different points in time (starting from 1874) through digital elevation models (DEMs), the latest of which refers to the year 2007. Vertical accuracy of the DEMs before and after 1960 has been assessed to be in the order of ±2 m and ±0.5 m, respectively (Bauder and others, 2007). The bedrock topography is known from both ground-based (Farinotti and others, 2009a,b) and helicopter-borne (unpublished data) radio-echo soundings. A mass-balance time series in daily resolution was reconstructed for the glacier starting from 1865 by Huss and others (2008), who combined direct stake measurements with geodetically derived ice volume changes. In the same study, representative time series of daily mean air temperature and daily precipitation sums were reconstructed for the glacier for the period 1865–2007. Farinotti and others (2012) extended the time series to 2010 using the same methodology. In the following, these time series are referred to as ‘observed meteorological time series’ as they rely entirely on direct observations retrievable from the archives of the Federal Office of Meteorology and Climatology (MeteoSwiss). The mean annual cycle of temperature and precipitation during the reference period 1980–2009 is shown in Figure 1c. According to the data, which refer to the elevation of the glacier center point (2700 m a.s.l.), mean annual temperature was −2.2°C, and mean annual precipitation 1940 mm.

Simulated meteorological data for the period 1961–2099 with daily resolution are available from the results of the ENSEMBLES project (Van der Linden and Mitchell, 2009). This European project, led by the UK Met Office, had the aim of providing a range of projections for future climate evolution including probabilistic information. For the analyses, the mid-range A1B emission scenario (Nakićenović and Swart, 2000) was used. Out of the ENSEMBLES database, six different ‘GCM/RCM model chains’ (i.e. combinations of global and regional climate models (GCMs and RCMs, respectively)) are considered Table 1). All model chains provide data at a horizontal resolution of ~25 km, and daily time series for the study area were derived by averaging the four model gridcells located closest to it. As described further below, these time series are only used to define plausible scenarios of altered variability, and not to force the ice-dynamics model directly, thus avoiding the eventual necessity of a bias correction.

Table 1. Overview of GCM/RCM model chains considered in this study. Acronyms for institutions: C4I, Community Climate Change Consortium for Ireland; DMI, Danish Meteorological Institute; ETHZ, Swiss Federal Institute of Technology, Zürich; METO-HC, Met Office Hadley Centre for Climate Prediction and Research, United Kingdom; MPIMET, Max Planck Institute for Meteorology, Germany; SMHI, Swedish Meteorological and Hydrological Institute

Methods

Statistical models for climate variability

There is no agreed definition of ‘variability’ for a particular climate variable. Here a framework which aims at decomposing the temporal fluctuations of a given meteorological time series into three components, called the year-to-year, month-to-month and day-to-day variability, is proposed. The decomposition permits these components to be acted on individually when using a weather generator to derive synthetic time series. The series are then used to force a mass-balance and ice-dynamics model, thus revealing the effect on the glacier evolution. The general idea is to consider deviations from a detrended long-term mean, describe these deviations with an appropriate statistical model, and parameterize the residuals of that model in order to obtain a suitable metric for variability. Detrending of the time series is necessary since ‘naively analyzing time series which have a trend leads to a trend-induced inflation of variability’ (Scherrer and others, 2005). Here a statistical model is defined as ‘appropriate’ if the residuals ε of that model have the properties of white Gaussian noise, i.e. are independent and identically distributed (i.i.d.) following a normal distribution with zero mean and scale parameter σ (i.e. i.i.d.). The scale parameter σ is then the desired ‘suitable metric for variability’.

For the sake of simplicity, and considering the information readily available in future projections, the analyses are constrained to the two variables temperature and precipitation. In the measured meteorological time series, the correlation between these two variables is weak at any time aggregation (Fig. 3). If any, a weak anticorrelation (coefficient of determination ~0.11) can be found for yearly values (very dry years tend to be warm; very wet years tend to be cold). This allows the two variables to be considered separately in the first instance, although some dependency will be introduced.

Temperature variability

Detrending of the temperature time series is performed with a STL decomposition (i.e. Seasonal Trend decomposition procedure based on Loess; see Cleveland and others, 1990). The time series detrended in this way is then aggregated to yearly values T _yr. Analysis of the autocorrelation function (ACF) and the partial autocorrelation function (PACF) of this new time series reveals that an autoregressive moving-average (ARMA) model (Whittle, 1951; Box and others, 2008) of order p = 1 and q = 1 is appropriate for describing it. In general, an ARMA(p, q) model describes the evolution of a given random variable X as

(1)

where X_t is the value at time t, c is a constant (c = 0 throughout this study), a_i and b_j are coefficients to be estimated and i.i.d. is a noise term with value ε_t at time t (the residual). The coefficients a_i and b_j are estimated according to the algorithm by Gardner and others (1980) as implemented in the software package R, and the values are listed in Table 2. For convenience, the ‘ARMA(1,1) model fitted to the time series of annual anomalies of the detrended time series’ will be referred to as the model. The scale parameter σ of the residuals of this model is estimated with (where sd(ε) denotes the standard deviation of ε) and is chosen to define the year-to-year variability. In the following, this parameter is denoted by (Table 2).

Table 2. Coefficients (coeff) and standard errors (s.e.) as estimated for the different ARMA models from the analysis of the observed meteorological time series during the reference period 1980–2009

The mean monthly temperature for month i (T _{mon, i} ) can be described in terms of a deviation ΔT _mon, _i from the mean annual temperature T _yr (i.e.T _{mon, i} = T _yr +ΔT _mon, _i ). This deviation can be further decomposed into a month-specific mean deviation which characterizes the mean deviation of month i from T _yr (e.g. on average, the mean temperature for January will be some degree lower than the annual mean, whereas, on average, August will have a mean temperature higher than the annual mean), and an additional deviation characterizing the deviation of the particular month from (e.g. the mean temperature of a particular January may well differ from the temperature expected considering the mean annual temperature for the given year and the average deviation of the month January from it). This concept, which is visualized in Figure 2a, can be formalized as

(2)

Fig. 2. Conceptual representation for the decomposition of the detrended (a) temperature and (b) precipitation time series. In each plot, a randomly selected time series of 1 year’s length is shown. Note that for precipitation, daily values are shown in the inset of the bottom right only. The individual symbols are explained in the text.

Fig. 3. Coefficient of determination for the dependency between temperature and precipitation in the observed (red) and simulated (black) time series during the reference period 1980–2009. For the simulated time series, the median (solid line) and the 95% confidence interval (dotted) are shown. Aggregations are performed as sum (precipitation) or mean (temperature) of daily values.

For the definition of month-to-month variability, a model for is required. Analysis of the ACF and PACF indicates that an ARMA(3,1) model is suitable for describing this term. The coefficients of this model, called the model, are estimated in the same way as for . Similarly, the scale parameter of the model residuals, , is the desired metric for the month-to-month variability (Table 2).

Exactly the same approach is used to define the day-to-day variability. In analogy to Eqn (2), one can write

(3)

Again, the concept is visualized in Figure 2a. For , inspection of the ACF and PACF indicates an ARMA(2,4) model. This model is called the model, and the scale parameter of the residuals, , is the metric for the day-to-day variability.

With the aim of mimicking some dependency between temperature and precipitation, the residuals of the model are analyzed separately for wet and dry days. The observation that the distributions differ in the two cases (not shown) is included by estimating two separate coefficients and for wet and dry days respectively (Table 2). This will be taken into account when generating synthetic time series. Note that by introducing the dependency at the daily timescale, the dependency is also propagated to any aggregated data.

Precipitation variability

Because precipitation cannot be negative, the framework as described for temperature cannot be applied in exactly the same way. However, the general idea is maintained. Since the annual precipitation sum P _yr can be regarded as strictly positive (years without precipitation can be virtually excluded in the considered climate setting), a log-transformation is suitable for deskewing the distribution (Mosteller and Tukey (1977) call this a ‘first-aid transformation’). After removal of the mean value, analyses of the ACF and PACF reveal that no particular model is required for having a signal with the characteristics of white Gaussian noise, and one can write i.i.d. The parameter for the year-to-year variability is then , and is estimated with (Table 2).

For the monthly values, a decomposition similar to that for temperature is applied, whereby the deviations are considered from the monthly sum which would arise if the yearly precipitation P _yr were uniformly distributed over the 365 days of the year. In analogy to Eqn (2) one can write

(4)

Here P _mon, _i is the precipitation sum for month i, n _days, _i the number of days in that month, the month-specific average deviation from the monthly sum given by the uniformly distributed precipitation, and the additional deviation for the particular month (Fig. 2b). Fitting a model directly to is not suitable, since the distribution of the quantity is slightly skewed (negative deviations are constrained by the fact that P _mon, _i ≥ 0). This problem is handled by adding the absolute value of the maximal possible negative deviation (in order to obtain positive values) and log-transforming the corrected values (for deskewing the distribution). The variable transformed in this way can be described with an ARMA(1,1) model (the model), and the scale parameter of the residuals, , is the metric describing the month-to-month variability. The coefficients estimated for the model are listed in Table 2.

Since precipitation is a positive quantity, the idea pursued so far is no longer suitable for describing daily values, and the approach is adapted. In particular, the two models set up for P _yr and P _mon are not required for the analysis of daily precipitation, but are used later when generating synthetic time series. The idea (first presented by Richardson (1981) and often used) is to describe the daily precipitation time series with a combination of a two-state first-order Markov chain (e.g. Deily, 1966) for the succession of wet and dry days, and a Γ-distribution for the precipitation sums in wet days. ‘Two-state first-order Markov chain’ means that the model (hereafter referred to as the Markov chain model, MCM) describes a process with only two possible states (wet or dry); the state for a particular time-step is only dependent on the state of the time-step before (‘first-order’), and the transition between states can be described through distinct probabilities (the idea of a ‘Markov chain’). Defining a state variable Z_t which assumes the values Z_t = 0 for a dry day and Z_t = 1 for a wet day, the (conditional) probabilities that need to be defined are (1) P00 = P{Z_t = 0 |Z_{t −1} = 0} (probability of having two consecutive dry days), and (2) P11 = P{Z_t = 1|Z_{t −1} = 1} (probability of having two consecutive wet days). The complementary probabilities are then given by P01 = P{Z_t = 1|Z_{t −1} = 0} = 1 – P00, and P10 = P{Z_t = 0 |Z_{t −1} = 1} = 1 – P11, and this fully describes the model. For the Γ-distribution, two parameters κ. and θ, describing the shape and scale respectively, need to be estimated. This is done with and where is the daily precipitation amounts for wet days. For the definition of day-to-day variability, P00 and P11 rather than κ. and θ are used, which seems the more reasonable choice since altering one of the two probability parameters necessarily influences the distribution of the daily precipitation amount as well. This is because when generating synthetic time series, the monthly precipitation sum will be given a priori by the model. The values estimated for the individual parameters are listed in Table 2.

Generation of daily weather data

The decompositions described above can be used in an inverse way to generate synthetic time series. In order to preserve the (weak) dependency between temperature and precipitation, precipitation time series are generated first. This is done in four steps: First, a time series with the desired length is generated for the yearly precipitation totals by taking a random sample from the distribution , and back-transforming the values (yearly sums were centered and log-transformed before). In a second step, a series of monthly anomalies of the given length is generated with the model. After back-transformation, these anomalies are superimposed on a year-by-year basis onto the monthly sum which is given by distributing the yearly precipitation uniformly over the days of the year. In rare cases where the superposition leads to negative values, the difference is uniformly distributed over the remaining month of the year (e.g. if for the October of a given year the superposition leads to a monthly precipitation sum of −11 mm, the precipitation sum for October will be set to 0 mm and the sum of all other months in that year decreased by 1 mm). This procedure guarantees consistency with the previously generated yearly totals in all cases. In a third step, a daily sequence of wet and dry days is generated with the MCM model, and in the fourth, final step, a daily precipitation sum is assigned to each day defined as wet. This happens by sampling a series of values from the fitted Γ-distribution. Consistency between daily totals and monthly sums is enforced by scaling the daily values with an appropriate factor.

Once a precipitation time series has been created, the associated temperature time series is generated according to the following scheme: First, a time series of yearly temperature anomalies is generated with the model. Second, the same is done for the monthly anomalies by using the model. In a third step, two time series are generated for the daily anomalies with the model. In the first time series, the variability parameter is set to , and in the second to . Consistency between the two time series is guaranteed by setting the random seed of the two generating processes to the same value. In a fourth step, the different anomalies are superimposed on the corresponding long-term means, giving a synthetic, detrended time series in daily resolution. When superimposing the daily anomalies, the value is chosen from one or the other time series according to the wet/dry condition given from the previously constructed precipitation time series. Finally, a trend is added to the time series in order to shift the mean value to the desired level (following the previous steps only would give time series with a long-term mean of 0°C). Note that within this framework, ensembles of time series generated with different sets of variability parameters have exactly the same mean values on a yearly, monthly and (for temperature) even daily basis.

Validation of the weather generator

The weather generator is validated by generating a sample of 100 time series for the reference period 1980–2009 and comparing a set of characteristics to the values derived from the observed time series. The set of characteristics includes: (1) For temperature: average and standard deviation (SD) of daily temperatures in the course of the year (365 comparisons for each time series, one per day of the year); average and SD of monthly temperatures in the course of the year (12 comparisons for each time series, one per month of the year); average and SD of annual temperatures (1 comparison per time series); average and SD of positive degree-days (PDDs) per year (i.e. yearly sum of temperatures >0°C); average warm-spell length (i.e. number of consecutive days with temperatures >0°C); and average cold-spell length (i.e. average number of consecutive days with temperatures ≤−5°C). (2) For precipitation: number of wet days per year; average and SD of precipitation amount for wet days; average and SD of monthly precipitation sums in the course of the year (12 comparisons per time series); average and SD of annual precipitation sum (1 comparison per time series); average and SD of the yearly precipitation occurring at a temperature <1.5°C (as proxy for accumulation); average wet-spell length (i.e. average number of consecutive days with precipitation); and average dry-spell length (i.e. average number of consecutive days without precipitation).

Following Kou and others (2007) and Min and others (2011), the performance of the weather generator is assessed by evaluating the number of individual, synthetically generated time series for which a significant difference occurs in a particular characteristic. The significance of the differences for mean and variances is assessed with the tests by Wilcoxon (1945) and Levene (1960), respectively. Spell lengths, PDDs and precipitation characteristics are log-transformed prior to testing for reducing the skewness. Differences are ascertained to be significant according to the P-value yielded by the test statistic, and the results grouped in four different classes of significance (Table 3) following Kou and others (2007).

Table 3. Validation of the weather generator. For each characteristic, the percentage of realizations generated for the reference period 1980–2009, and falling within a given difference level, is listed. Differences are classified to be very significant (VSD) for P < 0.01, moderately significant (MSD) for 0.01 ≤ P < 0:05, slightly significant (SSD) for 0:05 ≤ P < 0.10 and not significant (NSD) for P ≥ 0.1 (P being the P-value of the test statistic). The mean of the characteristic for the observed and synthetically generated time series is given by ‘obs.’ and ‘gen.’, respectively. Abbreviations: avg = average, sd = standard deviation, SL = spell length, WDY = number of wet days per year, PBT = yearly precipitation occurring at temperatures below the snow/rain threshold (i.e. 1.5°C). Temperature and precipitation statistics are given in 8C and mm, respectively, SL in days

The validation shows that the weather generator is capable of satisfactorily (i.e. ≥90% of the time series showing non-significant differences (NSDs) at the 10% level) mimicking the observed time series in most of the considered characteristics, and in the characteristics which are expected to have a major influence on the glacier response (e.g. the PDD statistics and the proxy characteristics for accumulation) in particular. Slight tendencies can, however, be observed in (1) underestimating the SD of daily temperatures if considered in the course of the year on a day-by-day basis (only 85.5% of NSDs), (2) underestimating the length of warm and cold spells (86% and 88% of NSDs, respectively), (3) underestimating the number of wet days per year (89% of NSDs), (4) overestimating the SD of the precipitation amount for wet days (78% of NSDs), (5) underestimating the SD of monthly precipitation sums if considered on a month-bymonth basis (87.1% of NSDs) and (6) underestimating the length of wet spells (87% of NSDs, which is a consequence of (3)).

One may argue that these tendencies can be regarded as acceptable for the presented application since (1) the differences occurring for the SD of daily temperature and precipitation are mainly contained in the 5–10% level of confidence, and (2) the role of precipitation will be shown to be negligible when compared to temperature, which relativizes the importance of the precipitation characteristics.

The degree to which the dependency between temperature and precipitation is preserved in the synthetically generated time series is assessed by computing the coefficient of determination (e.g. Box and others, 2008) between the two variables for aggregation lengths ranging from 1 to 365 days. Aggregation is performed by computing mean and total values for temperature and precipitation respectively. The obtained values are then compared to those calculated from the observed time series. Figure 3 shows that the characteristic of the dependency is well maintained. At most, the synthetic time series show a slight tendency to underestimate the dependency for long periods of aggregation.

Variability scenarios

For all models included in the weather generator, reference values for the variability parameters (i.e. , , P00 and P ₁₁, with i = [day, month, year]) are determined by analyzing the measured meteorological time series (last row in Table 4; red lines in Fig. 4). Scenarios for an altered climatic variability are then constructed in two different ways: In a first set of scenarios, the reference values are systematically doubled or halved (scenarios with ending ‘d’ in Table 4), and in a second set, the parameters are adjusted in order to explore the range of variability expected in the future according to the results of the ENSEMBLES model chains (scenarios with ending ‘E’ in Table 4). For the second set, the relative changes with respect to the reference period 1980– 2009 are computed for a given parameter and a given GCM/ RCM model chain in a running window of 30 year width (Fig. 4). The scenarios are then constructed in order to explore 95% of the range of these changes (bars on the right-hand side of Figure 4a–c and e–g).

Table 4. Summary of the scenarios (‘Scen.’) defined for temperature and precipitation variability. The σ-parameters are given as factors relative to the reference value determined from observed meteorological time series (‘Obs.’) while for the P-parameters the stated coefficients are additive. ‘Same as observed’ means that the time series of the variable was not generated with the given model but obtained by shifting the observed time series as described in the text

Fig. 4. Temporal evolution of the defined variability parameters (a–c, e–g), and long-term temperature and precipitation trends (d, h) according to six GCM/RCM model chains. Model chains with noticeable trajectories are labeled (see Table 1 for numbering). (a–c), (e) and (f) show relative deviations during a given 30 year interval from the values calculated for the reference period (1980–2009). In (d), (g) and (h) deviations are given as absolute values. The red curve depicts the values derived from the measured meteorological data. The bar on the right-hand side of (a–c) and (e–g) shows the range containing 95% of the data. In (d) and (h) the thick black line is a 30 year running median of the ensemble.

When generating the scenarios, three cases are distinguished depending on how the temperature and precipitation time series are obtained. In the first case, temperature is generated for the period 2010–2100 (91 years) with the corresponding models, while the time series for precipitation is constructed by adopting the observed records for the period 1920–2010 (also 91 years), adjusted with the long-term trend given by the GCM/RCM model chains (a more detailed explanation follows shortly). In the following, this set of scenarios will be referred to as the ‘temperature-only scenarios’ since the variability in precipitation is left unchanged. In the second case, the observed temperature for the period 1920–2010 is corrected with the long-term trend given by the GCM/RCM model chains, whereas the precipitation time series is generated with the given models. This set will be referred to as the ‘precipitation-only scenarios’ since the temperature variability is unaltered. In the third case, both temperature and precipitation are generated synthetically, and the set is referred to as the ‘combined scenarios’.

The temperature-only scenarios include nine different combinations in which the individual variability parameters are varied systematically: scenario T.ref is the reference scenario in which all parameters concerning the temperature time series are kept at the same level as estimated from the observed meteorological time series, scenarios T01–T06 are constructed by varying one variability parameter at a time, and in scenarios T07 and T08 all parameters are varied together. The 11 precipitation-only scenarios are constructed in a similar way. P.ref is the reference scenario, P01–P08 are scenarios in which the variability parameters for precipitation are varied one by one, and P09 and P10 are scenarios in which all parameters are varied together. The combined scenarios include scenario TP1, in which all variability parameters of both temperature and precipitation are increased, and scenario TP2, in which all variability parameters are decreased. The values for the individual variability parameters used in the different scenarios are reported in Table 4. Note that a particular case is given by the parameters steering the day-to-day variability of precipitation, since they express a probability. In this case, the range in which the parameters are varied in the scenarios is not defined through the relative deviation from the reference period, but by the spread of the parameter values in the six model chains. The spread is computed, and then centered around the parameter value estimated from the observed meteorological time series.

Since all scenarios are meant to represent possible future conditions, long-term trends are added to both the temperature and precipitation time series. The trends are defined to be the median of the deviations from the reference period computed for the six model chains over a running window of 30 years (black line in Fig. 4d and h). Linearly extrapolated to the period 2085–2114 and compared to the reference period 1980–2009, the trends result in an increase of +4.2°C in mean temperature, and a decrease of 51 mm in mean annual precipitation. These trends are added to all the constructed scenarios on a year-to-year basis, meaning that the time series of a particular year is generated first, and the value corrected a posteriori with the according trend. For temperature, the correction is additive (e.g. for the year 2100, 4.2°C are added to the generated time series), whereas for precipitation the correction is multiplicative (e.g. for the year 2100, the generated time series is multiplied by (2260 – 51)/2260b = 0.97, since the annual precipitation in the reference period was 2260 mm). For each of the named scenarios, a set of 100 synthetic time series spanning the time period 2010–2100 is generated, and used to force the combined glacier mass-balance and ice-dynamics model described hereafter.

Modeling of glacier evolution

The evolution of Rhonegletscher is simulated using the generated meteorological time series to force the ice-dynamics model presented by Jouvet and others (2008). The latter works in combination with a mass-balance model based on the approach by Hock (1999). This model combination has been used before to successfully reconstruct and project the transient evolution of both Rhonegletscher (Jouvet and others, 2009) and Grosser Aletschgletscher, Swiss Alps (Jouvet and others, 2011b), during the periods 1874–2100 and 1880–2100, respectively.

As described in Jouvet and others (2009), the ice-dynamics model treats the glacier ice as an incompressible non-Newtonian fluid and neglects inertial terms. Thus, the mass and momentum equations reduce to a stationary nonlinear Stokes problem in the ice domain Ω:

(5)

where ε(u) = (∇u + ∇u ^T )/2 is the rate of strain rank 2 tensor, ρ = 900 kg m⁻³ is the ice density and g = (0, 0, 9.81) m s⁻² is the vector of gravitational acceleration. Ice viscosity μ is a function of the velocity field u = (u, v, w), as described by Glen’s flow law (Glen, 1955):

(6)

in which n = 3 is the flow law exponent, A is the flow rate factor and bar is a regularization parameter which prevents infinite viscosity for zero strain. Boundary conditions for Eqn (5) are given at the glacier surface and base. At the surface, no forces apply:

(7)

(n is the unit outer normal vector along the boundary of the ice domain Ω, and p is pressure), whereas at the base the sliding condition reads:

(8)

In Eqn (8), t ₁ and t ₂ are two orthogonal vectors tangent to the ice–bedrock interface, and α a sliding coefficient given by

(9)

where C is a positive value to be tuned, and s ₀ = 0:01 m a⁻¹ is a regularization parameter which prevents infinite α for zero velocity.

Glacier evolution is then calculated by solving the transport equation

(10)

where φ = φ(x, y, z, t) is the characteristic function (φ = 1 if location (x, y, z) ∊ Ω at time t, and φ = 0 otherwise), and bρ is a source term describing the change of glacier mass due to accumulation and ablation (b is mass balance).

For any location (x, y, z), daily mass balance b is computed by subtracting daily surface ablation c, calculated as

(11)

(Hock, 1999), from daily surface accumulation a, calculated as

(12)

(Huss and others, 2008; Farinotti and others, 2012), i.e. b = a − c. Contributions of internal and basal mass balance are neglected. In Eqn (11), f _M is a melt factor,r _snow/ice are two distinct radiation factors for snow and ice, I _{pot, (x, y, z} ) is the potential direct clear-sky solar radiation at location (x, y, z) and is the mean daily air temperature at the same location. For days with , no ablation occurs. The spatial distribution of is obtained by interpolating the temperature synthetically generated for the reference location using a constant lapse rate. In Eqn (12), P _ref is the precipitation given by the synthetically generated time series (which refers to the center of the glacier, as a reference point ‘ref’), C _prec is a correction factor accounting for the gauge-catch deficit (Bruce and Clark,1966), z – z _ref is the elevation difference between the reference and the considered location, dP/dz is the lapse rate with which precipitation increases with elevation (Peck and Brown, 1962), D _{snow, (x, y)} is a spatially distributed factor which accounts for snow redistribution processes (e.g. Tarboton and others, 1995; Huss and others, 2009) and r _s is the fraction of solid precipitation. D _{snow, (x, y)} is determined from characteristics of the surface topography (Huss and others, 2008), while r _s decreases linearly from 1 to 0 in the temperature range T _thr − 1°C to T _thr + 1°C, where T _thr = 1.5°C is a threshold temperature distinguishing snow from rainfall (Hock, 1999).

Numerical solutions to the equations concerning the ice-dynamics model are found by decoupling the computation of φ from that of u and p through a time discretization, and using two different space discretizations to solve the transport problem (Eqn (10)) and the nonlinear Stokes problem (Eqn (5)). The transport problem is solved using a structured grid of small cubic cells (with the goal of minimizing numerical diffusion), whereas the nonlinear Stokes problem is solved on an unstructured mesh of tetrahedrons with larger size (since the task is computationally expensive). Further details on the numerical implementation of the model are provided by Jouvet and others (2008).

For Rhonegletscher, both the parameters of the ice-dynamics model and the parameters of the mass-balance model have been calibrated before by Jouvet and others (2009) and Farinotti and others (2012), respectively. For the present study, these values are adopted without modification (Table 5). Details of the calibration procedures, which include comparison to measured surface velocity fields, surface mass balance, ice volume changes and runoff, are provided in the aforementioned publications.

Table 5. Parameters (‘Param’) for the ice-dynamics and mass-balance model calibrated for Rhonegletscher (Jouvet and others, 2009; Farinotti and others, 2012). The column ‘Eqn’ indicates the equation in which the parameter appears

Results

The glacier evolutions resulting from the different model runs are shown separately for the sets of scenarios of doubled/halved variability (Fig. 5) and the scenarios based on the ENSEMBLES model chains (Fig. 6). For each individual scenario, the evolution of the total glacier ice volume (Figs 5a and b and Fig6a and b), a longitudinal profile of the glacier by the year 2055 (Figs 5c and d and 6c and d), as well as box plots for the distribution of the total ice volume (Figs 5e and g and 6e and g) and for the position of the glacier terminus at the same time (Figs 5f and h and 6f and h), are shown.

Fig. 5. (a, b) Temporal evolution of the total ice volume according to the different scenarios of doubled/halved variability (gray). At each point in time, results are given as the median of the 100 realizations comprising the scenario. The reference scenarios (scenarios T.refand P.ref for (a) and (b) respectively) are shown in red. Scenarios that noticeably depart from the reference are labeled. (c, d) Longitudinal glacier profile (see Fig.1a for location) as modeled for the year 2055. The glacier surface (gray) at any location is given as the median of the 100 realizations the scenario is composed of. The reference scenarios (red) and the glacier bedrock (black) are shown. (e, g) Box plots of the total ice volume calculated for the year 2055 according to the different scenarios. (f, h) Box plots of the glacier terminus position for the individual scenarios as calculated for the year 2055 (position is given as distance from the location of the 2010 glacier terminus). Box plots show minimum and maximum values (crosses), 0.025 and 0.975 quantiles (whiskers), 0.25 and 0.75 quantiles (box), and average (line inside box). Temperature-only and combined scenarios are shown in the left set of panels, precipitation-only scenarios in the right set. Table 4 gives details of the scenarios.

Fig. 6. Same as Figure 5, but for scenarios derived from the analysis of the ENSEMBLES model chains. Scale of the individual panels is kept the same as in Figure 5 for better comparison.

Three observations immediately catch the attention: (1) Large differences occur for the temperature-only and combined scenarios when the variability is doubled/halved (left panels in Fig. 5), (2) almost no effect is visible in the precipitation-only scenarios (right panels in Figs 5 and 6) and (3) no significant differences are visible for the scenarios derived from the ENSEMBLES model chains (Fig. 6).

The differences in the evolution of glacier volume for a doubled/halved temperature variability are especially striking: after 45 years of simulation, the total glacier volume for the T07d scenario (doubled variability at all three time-scales) is 63.8% smaller than in the reference scenario T.ref (Fig. 5a), whereas it is 18.4% larger for the T08d scenario (halved variability at all timescales). The effect is slightly less pronounced for the combined scenarios TP1d and TP2d (simultaneous doubling/halving in both temperature and precipitation variability), where the average difference in total volume is −61.6% and 14.6%, respectively. At this stage it may be worth recalling that for each scenario, the ensemble of meteorological time series used to force the ice-dynamics model has the same mean value on a yearly, monthly and (for temperature) even daily basis.

It is interesting to note that the effect of an altered temperature variability increases as the timescale of variability changes decreases. In fact, altering the year-to-year variability has less effect than altering the month-to-month variability (cf. results for scenarios T01d and T02d, or T04d and T05d), and similarly, altering the month-to-month variability has less effect than varying the day-to-day variability (cf. results of T02d and T03d, or T05d and T06d). On the other hand, the effect of a simultaneous variation of the parameters is, in first approximation, additive (cf. T07d and T01–03d, or T08d and T04–06d).

Differences in terminus position are also important (Fig. 5c and f), although the particular bedrock geometry, showing several prominent overdeepenings (e.g. around km3 and 4 in Fig. 5c), keeps the glacier tongue relatively stable at these particular positions. After 45 years of simulation and compared to the reference scenario T.ref, the terminus position for scenarios T07d and TP1d differs on average by ~0.5 km from the reference scenario, which is significant at the 5% confidence level. Concerning the particular bedrock geometry, it is worth mentioning that the configuration is likely to cause the formation of proglacial lakes, as is observed for the current (2013) terminus position. This would cause the glacier front to calve, and thus to retreat faster than predicted by the model. Such effects are, however, not taken into account in the current model setting.

Compared to temperature, changes in precipitation variability have only a minor effect. The largest difference compared to the reference scenario P.ref is found for the scenario in which only the year-to-year variability is varied (P01d). On average, the effect on the computed total volume by the year 2055 is in the order of 8.9% but the difference is not significant due to the large spread (Fig. 5b and g). The differences in terminus position are not significant for all the considered scenarios (Fig. 5h). The most noticeable feature is the increase in variance for the scenarios in which the year-to-year variability is increased (scenarios P01d and TP1d).

Very minor differences occur between individual scenarios derived from the analysis of the ENSEMBLES model chains (Fig. 6). This is true for both the temperature- and precipitation-only scenarios. Largest differences are found for scenarios T07E (increased temperature variability at all timescales) and T08E (decreased temperature variability at all timescales). In these cases, the total volume after 45 years of simulation differs on average by −7.7% and +5.1 %, respectively, from the reference scenario T.ref. In the evolution of the terminus position, the differences are even smaller: the average terminus position in scenario T07E differs from the reference scenario by only 0.1 km, which is not significant. Again, the bedrock geometry plays a role in dampening the differences.

Discussion

Especially for the temperature-only scenarios, the result that changing the day-to-day variability has a larger effect than changing the variability on longer timescales may look surprising at first glance. Similarly, it may not be immediately clear why changes in temperature variability have an effect while changes in precipitation variability do not. The explanations are, however, remarkably simple.

For the temperature-only scenarios, an explanation can be found by considering the cumulative PDDs (as proxy for the energy available for melt) to which the glacier is subject during the modeling period (2010–2100). Figure 7a shows how the cumulative PDDs are affected most by variations in the day-to-day variability. This is because with the proposed definitions an increased day-to-day variability causes the largest absolute temperature fluctuations (cf. observed , and in Table 4), which in turn causes the temperature to rise above the melting point more often than in the case of an altered year-to-year variability. A rather sloppy yet simple formulation of the above consideration is that ‘for glacier mass balance, a rather cold day followed by a rather warm day has a very different effect from a very cold day followed by a very warm day, although the mean temperature may well be the same’, and this is especially true when the fluctuations occur around the melting point. The incorporation of variance in PDD models goes back to Braithwaite (1985), who proposed a probabilistic approach to account for short-term temperature fluctuation in a model driven by monthly mean values, but the effect of an altered climatic variability has not been quantified so far, and may be surprising in its amplitude. Recall, however, that the mass-balance model used in this study is based on a degree-day approach, thus neglecting energy fluxes. In particular, neglecting the heat fluxes between the glacier surface and the first subsurface layers may lead to an overestimate of the effect of a higher temperature variability, since the cooling effect of very cold days is not considered.

Fig. 7. (a, c) Box plots of the cumulated PDDs to which the glacier is subject in the individual scenarios during the simulation period 2010–2100. All precipitation-only scenarios (PXX) are constructed with the same temperature time series. (b, d) Box plots of the cumulative June–August precipitation sum (P(JJA)) occurring at a temperature T < 1:5°C for the different scenarios and the simulation period 2010–2100. Table 4 provides details of the scenarios.

For precipitation, the answer is given by the fact that accumulation is a cumulative process, and that the timing of the precipitation occurrence does not play a role, as long as the phase of precipitation (meaning the repartition between snow and rain) remains unchanged. In this case, a sloppy formulation that clarifies the sentence could be ‘it doesn’t matter if it snows a bit today and a bit tomorrow, or nothing today and a lot tomorrow; accumulation will remain the same’. Within the presented framework, a change in precipitation variability does not change the total precipitation sum by definition and thus, as long as the temperature variability is kept at the same level, no changes in the repartition of the precipitation phase will be observed. This is because an altered precipitation variability will, on average, not shift the timing of precipitation occurrence, whereas an effect could be visible if, for example, changing the variability caused more precipitation to take place during the summer months. Similarly to temperature, however, it must be noted that the use of a degree-day approach neglects the components of the energy balance. In particular, the model does not account explicitly for changes in incoming solar radiation that are likely associated with a change in precipitation variability and, thus, cloudiness. At most, such changes can be considered to have been taken into account implicitly through the correlation between incoming solar radiation and temperature (cf. Eqn (11) and Ohmura, 2001), since the proposed framework allows a feedback between precipitation occurrence and temperature distribution (see ‘Generation of daily weather data’).

The question remains about the interplay between simultaneous changes in the variability of both temperature and precipitation. One may potentially argue that an increase in temperature variability may lead to an increased frequency of summer-snowfall events (SSEs) because of an increase in the frequency of cold spells, and that this could have an effect through the albedo feedback of the freshly fallen snow. Such an effect is indeed visible, but has only a minor influence on the glacier evolution according to the performed simulations (scenarios TP1 and TP2 in Figs 5 and 6).

Looking at the cumulative (i.e. added over the modeling period) precipitation occurring during the months June–August at a temperature below 1.5°C (as proxy for SSEs), it becomes evident that the occurrence of SSEs is controlled by the variability in temperature rather than in precipitation (Fig. 7b). The largest effects are given by the scenarios T07 (increased temperature variability at all three timescales, unaltered precipitation variability) and TP1 (decreased variability of both temperature and precipitation at all timescales). It is interesting to note that in the case of increased temperature variability, which increases the occurrence of SSEs because of the increased probability of cold summer days, increased precipitation variability has an attenuating effect (i.e. the difference of SSEs compared to the reference scenario is decreased), while in the case of a reduced temperature variability, which decreases the occurrence of SSEs, a reduction in precipitation variability increases the effect (i.e. the difference to the reference scenario is amplified). This is because increased precipitation variability makes the distribution of daily precipitation amounts more uniform (since the monthly totals are fixed for the mean of the ensemble), thus decreasing the likelihood of heavy SSEs (which occur more often when the monthly precipitation is clustered in a small number of individual precipitation events). The occurrence of SSEs is, however, reflected only weakly in the glacier evolution, since the accumulation given by SSEs during the whole modeling period is, in all scenarios, <3% of the total accumulation only.

The presented analyses deliberately focused on variability changes at relatively short timescales (yearly to daily), thus neglecting variations associated with rearrangements of the atmospheric flow patterns that take place over longer periods. For the precipitation totals, for example, the North Atlantic Oscillation (NAO; e.g. Wanner and others, 2001; Hurrell and others, 2003, and references therein) has been shown to exert a considerable influence on the weather regime of western and northern Europe (e.g. Auer and others, 2001; Casty and others, 2005), and may thus require consideration when aiming at providing long-term projections. However, for the time being, the availability of corresponding scenarios is very limited, so they were discarded in this study.

Similarly, it must be stressed that the analyses focused on a particular aspect of the uncertainty linked to simulations of glacier evolution only. In this work, all factors besides those controlling the climatic variability (e.g. the parameters of the mass-balance and ice-dynamics model, the long-term trend in the climate time series, or the emission scenarios the climate projections are based upon) were assumed to be known and left unaltered. This does not mean that these factors can be considered free of uncertainty when carrying out projections; rather, keeping them constant allowed, in this case, isolation of the effect the study was focusing on.

Concerning the proposed framework for generating synthetic weather data, three further remarks are appropriate: (1) For precipitation, the fact that modeling the sequence of wet and dry days using two-state transition probabilities and deriving the daily rainfall amounts from a Γ-distribution is prone to underestimate the variance of daily precipitation has been noticed before (e.g. Srikanthan and McMahon, 2001). Although the implications were not serious in the present case (since precipitation variability was shown to have a minor effect on glacier evolution), more sophisticated methods would be required in further applications to better account for the clustered nature of the storm arrival process (Kilsby and others, 2007). (2) The proposed framework allows separate control of the variability at the yearly, monthly and daily timescales. The choice of these timescales may be criticized as somewhat arbitrary, since a breakdown into other aggregations (e.g. weekly instead of monthly, seasonal instead of yearly timescale) could, equally, have been considered. This criticism can certainly not be ruled out completely, but it seems hard to argue that another disaggregation would have been more justified. (3) The developed framework implicitly assumes independence of the variability at the three considered timescales. This design was intentionally chosen given the goal of the study (i.e. to explore the effect of variability at different timescales) but may be considered as a drawback in other applications.

The presented variability scenarios had a twofold objective: the first set of scenarios, in which the variability of the observed meteorological time series was doubled/halved, was included to show that changes in variability can potentially have very large effects in glacier evolution, whereas the second set of scenarios, obtained from the ENSEMBLES model chains, was included to answer the question whether such changes are likely to play a major role in the near future.

In this context, it is worth mentioning that recent RCM experiments agree in predicting increased variability during the coming decades. In a modeling study, for example, Schär and others (2004) found that European summer surface temperature variability might even double within the current century, while Scherrer and others (2005) stated that ‘although there are substantial differences between models regarding the location and amplitude of this effect, these results have qualitatively been confirmed by several other … studies (Meehl and Tebaldi, 2004; Brabson and others, 2005; Giorgi and Bi, 2005; Vidale and others, 2007)’.

The analyses presented here show that these changes are unlikely to have a significant effect on glacier evolution in the next century (the scenario derived from the ENSEMBLES model chains was constructed upon the 0.975 and 0.025 quantiles of the projected changes, and can thus be considered rather extreme), but suggest that at least caution is required when considering longer timescales. Reconstructions that span several centuries or millennia should at least discuss the possibility that the observed changes may also be linked to changes in climatic variability, rather than to the long-term trend of climate alone. This is true for reconstructions based on paleo-glacial landforms, which are directly linked to glacier evolution as addressed in this study, but may be similarly relevant for other analyses based on ice coring, dendrochronology or sedimentology, where the interpreted signals can show similar nonlinear effects to glacier mass balance.

Conclusions

Based on simple time-series analysis and on the decomposition of detrended temperature and precipitation series, a framework allowing the distinct definition of year-to-year, month-to-month and day-to-day variability has been presented. By forcing a three-dimensional full-Stokes ice-dynamics model set up for Rhonegletscher with ensembles of synthetically generated time series of temperature and precipitation, it was shown that glacier evolution can potentially differ considerably (up to 64% difference in ice volume within half a decade) depending on the variability scenario, but that for variability changes expected in the near future, the effect is likely to play a minor role (<8% even for extreme scenarios). This latter conclusion was drawn from a set of scenarios based on the output of six different GCM/RCM-model chains run within the ENSEMBLES project. The observation that different ensembles with the same mean temperature and precipitation evolution can lead to large differences in glacier evolution is explained by considering the sum of PDDs the glacier is exposed to. The effect of changed precipitation variability was shown to be of minor relevance. Although the interpretation of the results seems obvious once the explanations are given, the study shows that one should at least be cautious about making the common assumption of unchanged variability.