3.1. Bivariate formulation of Lliboutry’s model

To generalize Lliboutry’s model to the actually observed pair , we introduce spatial terms and temporal terms for winter balance under the additional constraint :

(2)

where superscript w denotes the winter period. The overall mean is the mean glacier-wide winter balance over the studied period, whereas the series represents the annual deviations to the mean. N ₂ denotes the two-dimensional (2-D) Gaussian vector of variance–covariance matrix ∑ which includes the marginal variances and . The residual pair reflects the discrepancy between the observations and the linear model which is expected to be low if the model fits the data properly. It is simply given by:

(3)

The main advantage of the joint formulation is that the covariance is also available, providing easy access to the correlation coefficient between the residual pair :

(4)

This makes recourse to an empirical estimation procedure such as that used by Reference Rasmussen and AndreassenRasmussen and Andreassen (2005) unnecessary.

The pair of annual deviations is assumed to be sampled from a 2-D Gaussian vector of variance–covariance matrix ∑ _β . Its mean is set to zero for consistency with the imposed constraints:

(5)

The hidden quantities have a hybrid status, i.e. they behave as observations with regard to the overall parameters, but they are parameters with regard to the observations. Such a multilayer construction characterizes a hierarchical structure. Its advantage over a non-hierarchical approach is to explicitly model the similarity of the different annual terms through the different terms of the ∑ _β covariance matrix. The marginal variances and quantify the magnitude of interannual fluctuations of respectively annual and winter deviations. As for the residuals, the covariance between annual and winter balance is given by , which provides easy access to their correlation coefficient . It is definitely non-negligible, since mass balance physically depends on winter balance:

(6)

The summer balance is obtained by subtracting winter balance from mass balance:

(7)

where superscript s denotes summer balance terms. Because of the properties of Gaussian linear models, it is also normally distributed:

(8)

and its variance decomposition into a spatial term , an annual term verifying and a residual is readily available:

(9)

Keeping an arithmetic averaging over the glacier surface, the mean summer balance over the studied period is then simply . Moreover, marginal N variances for the summer balance residual and for the annual deviation to the mean summer balance can be computed from the covariance matrixes ∑ and ∑ _β , respectively. This implies that the error made while estimating the summer balance using Equation (9) is higher than if mass balance and winter balance were independent, but remains smaller than σ _a+σ _w for σ _s, and smaller than for :

(10)

(11)

(12)

(13)

Finally, the correlation coefficients between summer balance and the two other terms are not explicitly available in the model and must therefore be estimated empirically. Since all residuals and annual deviations are, by definition, centred random numbers, the correlation coefficients and between the residuals pairs and are given by:

(14)

(15)

whereas the correlation coefficients and between the annual deviations are given by:

(16)

(17)

Equation (2) separates different independent effects using variance decomposition. We introduce the ratios λ _a, λ _w and λ _s (Equations (18–20)) to quantify the specific contribution (% of variance) to the total variability of the annual terms for annual, winter and summer balances respectively. Note that the denominator is an approximation of the total variability neglecting the spatial variability of the α_i terms. This is reasonable for a small glacier such as Sarennes (Reference Thibert and VincentThibert and Vincent, 2009), but would not be appropriate for a wider glacier covering a large elevation range:

(18)

(19)

(20)

Figure 2a summarizes this first model, noted M₁, using a direct acyclic graph (Reference LauritzenLauritzen, 1996). On the one hand, the parameters , i ∈ [1, N] and hidden variables have to be estimated using the data. On the other hand, all the other unknown quantities such as the summer terms, correlation coefficients and ratios are explicit combinations of and , which somewhat restricts the dimension of the estimation problem. With regard to Reference Rasmussen and AndreassenRasmussen and Andreassen’s (2005) approach, our model offers a fair quantification of the uncertainty levels thanks to the explicit modelling of the covariances ∑ and ∑ _β . With respect to Lliboutry’s model, it is also characterized by its fully bivariate nature. Figure 2b indicates how point variables are related to glacier-wide variables.

Fig. 2. (a) Direct acyclic graph of model M₁. Arrows express conditional dependence, circled nodes represent stochastic hidden variables, diamonds represent the overall parameters and rectangles indicate observed values. (b) Links between point and glacier-wide variables illustrated for the annual term. The alternative to estimate the mean annual glacier-wide balance from a geodetic balance (photogrammetry), Φ, instead of the spatial mean 〈B ^a〉 is also indicated.

3.2. Extracting structured temporal trends

With model M₁, the annual terms are modelled as if they were exchangeable, which means that they are assumed to be sampled from the same distribution whatever the considered year. This makes it possible to quantify the interannual variability, for instance thanks to the temporal ratios of Equations (18–20). On the other hand, a possible structured time trend that affects the whole glacier in a similar way and may therefore be related to climate change cannot be inferred. To overcome this limitation, we relax in this subsection the assumption of exchangeability of the . Moreover, for better consistency with the physical processes involved, we prefer working with the seasonal pair rather than with the pair .

A time-series analysis of the can theoretically be achieved using various time-series models developed in other fields, such as econometrics (Reference Booth and SmithBooth and Smith, 1982), hydrology (Reference Rao and TirtotjondroRao and Titrojondro, 1996), meteorology (Reference BloomfieldBloomfield, 1992) or climatology (Reference Von Storch and Zwiersvon Storch and Zwiers, 2002). The choice of the time-series model must be made with regard to the investigated temporal pattern, for example, memory effects, monotonic trends or change points for the mean (Reference Mearns, Rosenzweig and GoldbergMearns and others, 1997), variance (Reference DiazDiaz, 1982) or extreme values (Reference Nogaj, Parey and Dacunha-CastelleNogaj and others, 2007).

However, the balance between the available information and the number of model unknowns must also be considered to make it possible to extract the temporal structure from the hidden variables without producing an uncertainty level that is too high for any practical conclusions. We therefore restrict our purpose to a relatively simple change-point model with limited parameters. This model is relevant to Glacier de Sarennes data (see section 4) and has already been applied to similar problems in hydrology (Reference Perreault, Bernier, Bobée and ParentPerreault and others, 2000a,Reference Perreault, Bernier, Bobée and Parentb).

We designate (τ _w, τ _s) as the years of possible change points separating two periods of different conditions for winter and summer balances respectively. Before and after the change points, both the (Equations (21) and (22)) and the (Equations (23) and (24)) are broken down into random noises and structured linear trends. The random noises, with respective variances _: and , model the unstructured interannual fluctuations. The notation , for example, indicates that the parameter can take different values before and after the change point, so as to allow τ _s and/or τ _w to also reflect a change point in variance. The linear trends, expressed as w ^*(t) = a _w:+b _w: t and s ^*(t) = a _s: + b _s:t, are designed to capture gradual changes affecting summer and/or winter mass-balance terms. Seasonal components of the balance being expressed in m w.e. a⁻¹, the trend slopes (b _w., b _s.) are expressed in m w.e.a⁻². Depending on the continuity of w ^*(t) and s ^*(t) around τ _w and τ _s, the change points can be brutal, with a clear separation of two budget regimes, or just a distinction between two different trends (b _w1, b _w2) and/or (b _s1, b _s2). This model, noted M₂, is therefore flexible enough to capture monotonic trends and various types of changes in mean and variance for both summer and winter balances:

(21)

(22)

(23)

(24)

The trend for annual mass balance, a ^*(t), is easily accessible due to the linearity of the model:

(25)

If τ _s = τ _w, there are two linear trends b _a: = b _w: + b _s: and intercepts a _a: = a _w: + a _s: in mass balance. However, as soon as τ _s ≠ τ _w, three possibly different linear trends and intercepts are obtained for respectively t ≤ min (τ _s, τ _w), t ∊ [min (τ _s, τ _w) + 1, max (τ _s, τ _w)] and t ≥ max (τ _s, τ _w) + 1. This makes it possible to model different combinations of changes in winter and summer budgets, leading to three possible regimes of change in mass balance.

A major difference between models M₁ and M₂ is that, with model M₂, the variance of the annual deviations is no longer constant, depending on the existence of change points in the data. Nevertheless, it is possible to use Equations (26–28), where E denotes the mathematical expectation, to globally evaluate the dispersion of the annual terms in a way similar to that used in model M₁:

(26)

(27)

(28)

These global variances can be used to compute global temporal ratios quantifying the specific contribution to the total variability of the annual terms all over the studied period using Equations (18–20). Similarly, the global coefficient correlations and can be evaluated using Equations (16) and (17) respectively. For , however, Equation (6) is no longer valid, and empirical estimation must be used. This shows that, with regard to model M₁, the bivariate formulation at the annual level is lost with model M₂:

(29)

Finally, to compare the respective contributions to the annual variability of the temporal structure and the random fluctuations, we define the global ratios μ _a, μ _w and μ _s (Equations (30–32)). They are able to indicate whether the hypothesized structured linear trends with change points explain a large part of the annual variability of annual, winter and summer balances respectively, in a way similar to that of the classical determination coefficient R ² in a non-hierarchical linear regression:

(30)

(31)

(32)

Obviously, the different proposed ratios as well as the correlation coefficients between the can also be computed separately before and after the change points τ _w and τ _s instead of globally over the entire observation period. This is, for instance, straightforward for λ _s, since it only implies computing the ratio with the two values and However, for mass balance, the three possible periods must be considered, which multiplies the number of indicators to be interpreted without necessarily adding much information to the analysis. This is the main reason why we mainly consider the global quantities over the T years of study.

Summing up, by comparison with the first model, model M₂ involves the additional parameters and hidden variables w ^*(t), s ^*(t) that have to be estimated from the data. All of them are related to the quantification of the temporal structure, which is assumed to be well modelled by two linear trends possibly separated by change points. The rest of the model is not modified, except for the correlation structure at the annual level, which is no longer constrained by ∑ _{β 12}.

3.3. Bayesian inference and associated MCMC sampling

The different models proposed in this paper have been inferred within a Bayesian framework (Reference BergerBerger, 1985; Reference Parent and BernierParent and Bernier, 2007). The principle is to compute the joint posterior probability distribution function (PDF) of all parameters, hidden variables and missing values. Poorly informative marginal priors were used to ‘let the data speak for themselves’. For instance, pseudo-uniform priors were chosen for the spatial terms, poorly informative inverse gamma and inverse Wishart priors for the variances and variance–covariance matrixes, and poorly informative normal priors for the trend parameters.

Markov chain Monte Carlo (MCMC) schemes have been used to sample the joint posterior PDF. These powerful simulation methods are detailed, for example, by Reference Gilks, Richardson and SpiegelhalterGilks and others (1996) and Reference BrooksBrooks (1998). Examples of applications in a related field are given by Reference Eckert, Parent, Naaim and RichardEckert and others (2008, Reference Eckert, Naaim and Parent2010b). In our case, the efficient Gibbs sampler (Reference Geman and GemanGeman and Geman, 1984) could be run easily because of the Gaussian nature of a variance analysis that makes the full conditional distributions of the different parameters and hidden variables accessible.

Convergence of the MCMC sequence was checked by graphical comparison and various tests on the distributions obtained with different chains starting at different points of the parameter space (Reference Brooks and GelmanBrooks and Gelman, 1998). When the ergodic state was clearly reached, a sample of the joint posterior distribution was stored. For each quantity of interest, a point estimate and the related uncertainty was obtained by considering the marginal posterior mean and the related 95% credibility interval (i.e. the Bayesian counterpart of the classical confidence interval). Figure 3 shows an example of the marginal posterior distributions obtained for six quantities of interest. All are nicely shaped and unimodal, which indicates that information contained in the prior has been updated with the data and that the point estimates are meaningful. Among the distributions presented, those corresponding to σ _a and are direct results of the MCMC procedure. The others were obtained by computing their values at each step of the sampling procedure.

Fig. 3. Marginal PDFs for a few parameters, model M₁. (a) Standard deviation for the annual balance residuals. (b) Correlation coefficient between annual and winter balance residuals. (c) Standard deviations for annual balance deviations. (d) Correlation coefficient between annual and winter balance deviations. (e) Glacier-wide mean balance over the period of record. (f) Ratio between variance of annual deviations and total annual variance.

4.1. Lliboutry’s bivariate model

Table 1 summarizes the estimates obtained by applying model M₁ to the Sarennes series. It appears that 〈B ^a〉, the glacier-wide mean annual balance computed over the whole period under the assumption of equally weighted sites, is strongly negative. The glacier has lost ∼1 m w.e. a⁻¹ because, in mean, summer balance has been much larger in absolute value (−2.66 m w.e. a⁻¹) than winter balance (1.69 m w.e. a⁻¹). These results do not differ significantly from those obtained with the altitude–area integration relationship ; Table 1). They are also clearly consistent with previous studies of the same glacier (Reference VallaValla, 1989; Reference Vincent and VallonVincent and Vallon, 1997; Reference Thibert, Blanc, Vincent and EckertThibert and others, 2008; Reference Thibert and VincentThibert and Vincent, 2009) and other glaciers in the French and Swiss Alps (Reference Müller and KappenbergerMüller and Kappenberger, 1991; Reference VincentVincent, 2002; Reference Vincent, Kappenberger, Valla, Bauder, Funk and Le MeurVincent and others, 2004; Reference Huss and BauderHuss and Bauder, 2009).

Table 1. Posterior characteristics of model M₁ applied on Glacier de Sarennes point mass-balance series for the 1949–2007 period. SD denotes standard deviation, 2.5% and 97.5% denote the lower and upper bound of the 95% credible interval. For comparison, mean annual glacier-wide balances obtained with the glaciological method using the area–altitude distribution averaging (AAD subscripts) are also indicated for the annual, winter and summer terms

Marginal variances (σ _a, σ _w, σ _s) of the residuals are equal to 0.3, 0.21 and 0.43 m w.e. a⁻¹ for annual, winter and summer balance respectively, which is consistent with the inequalities of Equation (12). Marginal variances of the annual terms , , are more than three times higher. Moreover, spatial variability is relatively small (the values obtained on the five stakes are not very different from each other for all three variables; Table 2), which is consistent with the small size of the glacier. Note also that estimation error is small for the three α_i series, i.e. not more than 0.044 m w.e.a⁻¹. All this indicates that the interannual variability between the is the major source of variability at the glacier scale, constraining the three ratios λ _a, λ _w and λ _s to be above 90%. It also a posteriori justifies neglecting the spatial variability while computing the ratios of Equations (18–20).

Table 2. Spatial α _i terms, model M₁ applied on Glacier de Sarennes point mass-balance series for the 1949–2007 period. Indices 1–5 correspond to the five stakes located in Figure 1. Units are m w.e. a⁻¹

Residuals account for nonlinear effects that are not systematic in space and/or time and cannot therefore be captured in the α and β terms. Their correlations are all significantly nonzero at the 95% credibility level, which indicates that when the linear model performs poorly for one of the three terms of the balance, it is also biased for the two others. The positive values for and are easily understandable since when a local excess of winter mass input or summer mass loss occurs at a given stake in a given year, it influences accordingly the annual mass balance at this location. A nonzero value for is less intuitive. A physical explanation for a nonzero correlation between winter and summer balance residuals could be local excesses due to avalanches or local excesses/deficits due to snowdrift. In case of an excess, the corresponding albedo feedback in the summer balance would then result in a positive correlation. Indeed, an expanded snow ablation period would decrease ablation due to higher albedo at the glacier surface and reduced shortwave radiation absorption. Therefore, an enhanced (less negative) summer balance would be expected. In our case, is negative, and harder to interpret, but anyway the correlation remains low (−0.23).

At the level of annual terms, things are simpler: a very strong positive correlation exists between annual and summer balance. Identically, annual and winter balances are well correlated ( being significantly positive). Finally, winter and summer balances are uncorrelated with . This suggests that the albedo feedback mentioned above is not effective at the overall glacier scale.

In the Bayesian framework, the joint posterior distribution of the annual deviations is obtained. It is summarized in Figure 4 using box plots. The glacier-wide annual terms are obtained by adding the interannual mean, i.e. 〈B ^a〉,〈B ^w〉 and 〈B ^s〉 respectively, to each value (Fig. 2b). For all three series, the interannual fluctuations are very strong, so that the different annual terms are largely significantly different from each other at the 95% credibility level. This reflects the variability under unsteady climatic forcing conditions of the accumulation and ablation processes governing the two seasonal components of the balance. Even if apparently random fluctuations are strong, glacier-wide balance (i.e. ), although remaining almost always negative, seems to increase during slightly more than the first half of the studied period, and then to decrease again, reaching strongly negative values at the end of the study period. On the other hand, the winter mass input decreases over the entire study period, but with a strong shift in the middle. Glacier-wide summer balance is nothing more than the difference between the two other terms, with values much lower after 1985 than before.

Fig. 4. Box plots of β_t terms (annual deviations), model M₁. (a) Annual balance, (b) winter balance, (c) summer balance. The posterior mean is shown in black, the interquartile range in grey and the 95% credibility interval by dashed lines. Extremely similar plots are obtained with models M₁–M₅. As shown in Figure 2b, the glacier-wide terms can be obtained by adding to each β_t term the interannual mean given in Table 1, i.e. 〈B ^a〉,〈B ^w〉 and 〈B ^s〉 , respectively.

Even if the behaviour of the different annual terms is not identical, there is evidence in the three time series that changes separating different structured trends have occurred around the middle of the study period. This confirms the results of Reference Vincent and VallonVincent and Vallon (1997), reporting a strong positive accumulation change (+17%) on Sarennes since 1977 and a large increase in summer ablation and ablation rates (50%) since 1981. The latter has been thoroughly analysed by Reference Vincent, Kappenberger, Valla, Bauder, Funk and Le MeurVincent and others (2004). However, in these approaches, a single stake was used, and the significance of these changes was not assessed using specific statistical tests. Furthermore, annual terms were considered as observed, without taking into account the uncertainty related to their estimation. This justifies testing the hierarchical and non-exchangeable model M₂ on the same dataset in section 4.2 to identify temporal break points in means/trends/variances of the three annual terms.

For the three β_t series, estimation error measured as the width of the 95% credible interval for each estimate remains nearly constant over the whole period. The only exceptions correspond to the years for which one or several measurements are missing (e.g. 1951, three missing values out of five, or 1956, four missing values out of five), and for which estimation error is then logically higher. This ability to fairly reflect the level of available information is for us an important advantage of our modelling approach with regard to empirical estimation procedures. Furthermore, the uncertainty around the annual estimates remains tenable even with two or three missing measurements out of five. This is another clear benefit of a bivariate hierarchical model. Indeed, correlations between the different variables (Equation (2)) and from one year to another (Equation (5)) are explicitly modelled, so that information is transferred from one variable/year to another, thus reducing all credible intervals, and especially those corresponding to years with one or several missing values.

Similarly, missing data are estimated with credibility levels depending on the number of measurements for the considered stake/year. Figure 5a thus presents the annual balance at stake 3. Measurements exist between 1958 and 2007. Between 1949 and 1957, missing values are estimated on the basis of correlations between variables, sites and years. Variance of the estimated missing value is lower for 1957 than for 1956, since one more measurement is available (Fig. 5b and c).

Fig. 5. Estimation of missing data in the annual mass-balance series at stake 3 (missing values from 1949 to 1957). (a) Completion of the record at stake 3, , and associated uncertainty for the missing values. (b, c) Posterior distribution of mass balance at stake 3 estimated in (b) 1956 and (c) 1957.

Finally, model adequacy with respect to data can be checked by comparing the residuals with the centred linear model with variance covariance E. Figure 6 shows that, for the point mass-balance series, model fit is relatively good, with no bias and a nearly Gaussian distribution of the residuals, which is not surprising since the two-way ANOVA model has been known since Lliboutry’s work to be well suited to mass-balance series. Model fit is also not that bad for the winter point balances, which has been less often reported. However, with a statistical test such as the Kolmogorov–Smirnov test, the null hypothesis that both residuals are sampled from centred normal distributions of respective marginal variances and is questionable: it is accepted at the 95% significance level for annual balance, but only with a p value of 0.2, and rejected at the same significance level for winter balance with a p value of 0.02, just below the rejection threshold.

Fig. 6. PDF of residuals for linear model M₁ in annual (a) and winter (b) point mass-balance series.

These results indicate a small but significant discrepancy between the model and reality, especially for winter balance. In the latter case, it is probably related to the lower stability from year to year of the spatial structure of winter balance because of processes such as avalanches or snowdrifts whose influence has already been mentioned while discussing . The development of a ski resort on a part of the glacier and snow grooming may also have affected nonlinearly the spatial distribution of winter snow deposition. Another reason for the relatively low p values obtained is related to the normality test which is somewhat severe here as the quantity of available data is relatively large. This is not the case in other glaciological studies applying Lliboutry’s model which generally rely on much shorter measurement series. Note in any case that we are primarily interested in the annual terms that capture the main part of data variability rather than in the residuals, making the obtained low p values sufficient with regard to the scope of the study.

4.2. Change-point model

Model M₂ shares many common features with model M₁, so that some of the results are identical when both models are applied to the Sarennes dataset. For instance, this is the case for most of the quantities related to the observation layer: overall means, residual variances and correlation coefficients, as well as spatial terms. Their values (Table 3) therefore do not need to be discussed again.

Table 3. Posterior characteristics, model M₂ applied on Glacier de Sarennes point mass-balance series for the 1949–2007 period.

However, as soon as the temporal aspect is considered, differences appear between models M₁ and M₂. The introduction in model M₂ of a non-exchangeable structure and the related absence of an explicit specification of the correlation structure between the triplets explains why the overall variances and correlation coefficients differ slightly with model M₂ (Table 3) from their previous values (Table 1).

More precisely, all marginal variances are reduced, and no longer satisfies Equation (11). This is logical, since the means of the are, with model M₂, functions of time rather than set to zero, allowing better fitting to the data. As a consequence, the temporal ratios λ _w and λ _s are lowered by ∼10% compared to their values with model M₁, but the interannual variability remains largely predominant, still representing >80% of the total variability. For the correlation coefficients, the main features are slightly enhanced with regard to model M₁: and are more strongly positive, whereas there is still no correlation between winter and summer balance components. All these differences are, however, insufficient to make the time series from model M₂ graphically distinguishable from those of model M₁ displayed in Figure 4.

The values obtained for μ _a, μ _w and μ _s show that the assumed linear time trends represent ∼20% of the inter-annual variability of annual and winter balances, and 27% of the interannual variability of summer balance (Table 3). The years of change found for winter and summer balances are significantly different. The posterior distribution of τ _w shows a well-marked main mode in 1976, indicating strong evidence that a change actually occurred around that year for winter balance. For summer mass balance, the main mode is less marked in 1982, and other plausible years range between 1981 and 1984 (Fig. 7). This suggests that, for Glacier de Sarennes, the winter mass-balance regime changed ∼6 years before the summer balance regime.

Fig. 7. (a, b) Posterior distributions of the year of change for model M₂ for summer (a) and winter (b) balances. (c, d) Posterior distributions for change points in additional models M₃ (c) and M₄ (d).

is only slightly lower than , so τ _w cannot be considered as a change point in variance for the winter balance. On the other hand, for summer balance, variance around the interannual trend is much higher than . τ _s is therefore also a change point in variance separating two periods with interannual variabilities differing from the overall variance of Equation (28).

The mean trends and the related uncertainties are shown in Figure 8. The box plots may look quite wide, but this is because the linear fits are made on hidden variables, thus taking into account the uncertainty related to the for the estimation of the time trends. This may prevent firm conclusions, but once again properly reflects the amount of information present in the data. Logically, the uncertainty level is higher at the beginning and at the end of the study period where information density is lower, as well as around the change points. This is especially true for a ^*(t) which is influenced by both change points (τ _w, τ _s).

Fig. 8. Box plots of the structured trends, model M₂: (a) annual balance, (b) winter balance and (c) summer balance. The posterior mean is shown in black, the interquartile range in grey and the 95% credibility interval by dashed lines.

Looking closer, the two trend parameters (b _s1, b _s2) are significantly nonzero at the 95% credibility level, positive before the change point and negative after. On the other hand, the two b _w. trend parameters and the four intercepts are not significantly nonzero at the 95% credibility level (Table 3). However, even if checking that zero belongs to a credible interval is in some ways the Bayesian counterpart of a classical significance test, it is not a test in the strict sense of the term, and the shape of the posterior distribution must be taken into account to properly interpret the results. For instance, the marginal posterior distributions of the (b _w1, b _w2) trend parameters are far from centred around zero, indicating that there is a high probability that they may be nonzero. The decreasing trend before and after the strong shift τ _w must therefore be noted, even if it is definitely less marked than the two trends in the summer mass-balance series. Hence, the high uncertainty around the most probable values must be linked to the difficulty of extracting information from hidden layers rather than to the absence in the series of any temporal structure.

For winter balance, the structured trend is strongly discontinuous around τ _w = 1976, indicating that the change is sudden. Between 1949 and 1976, the mean glacier-wide winter budget , but is not uniformly distributed since there is a decrease of −0.37 m w.e. a⁻¹ over this period (b _w1 = −0.013 m w.e. a⁻²). In 1977, winter mass input increases strongly and its mean between 1977 and 2007 is 1.88 m w.e. a⁻¹.Note though that, as observed before 1976, winter mass balance slowly decreases in this second period by 0.20 m w.e. a⁻¹until 2007 (b _w2 = −0.007 m w.e. a⁻²). Mean glacier-wide winter balance was therefore lower than the long-term mean 〈B ^w〉 = 1.69 m w.e. a⁻¹ before the change point, and higher than the long-term mean after the change point, with the years of lowest values being just before the change point in 1976, and the years with the highest values just after.

Regarding summer balance, the transition lasts ∼4–5 years around τ _s = 1982. During these 4–5 years, it is shifted down by about 1.11 m w.e. a⁻¹. Before 1982, mean glacier-wide summer balance is , that is 0.4 m w.e. a⁻¹ more than the long-term mean 〈B ^s〉. Furthermore, during this first period, summer glacier-wide balance increases by 0.62 m w.e. a⁻¹ (b _s1 = +0.024 m w.e. a ⁻²). After 1982, summer balances are significantly lower, , i.e. nearly 0.6 m w.e.a⁻¹ below the long-term mean. During this second period, there is also a decreasing trend from year to year (b _s2 = −0.015 m w.e. a⁻²), with summer glacier-wide balance more in deficit by 0.30 m w.e. a⁻¹ between τ _s and 2007.

For glacier-wide annual balance, three regimes are separated by the two change points, τ_w and τ_s:

1. 1. Before τ_w = 1976, Glacier de Sarennes loses mass with a mean annual wide balance of −0.80 m w.e. a⁻¹, but the trend goes from strongly negative balances to less negative balances from year to year (b _a = +0.011 m w.e. a⁻²), mainly due to less negative summer balances.

2. 2. Between 1976 and 1982, glacier-wide winter balances increase strongly while the summer regime has not yet changed. As a result, the mean mass balance is slightly above 0 m w.e. a⁻¹ and the glacier, very close to equilibrium, maintains its mass. The trend is to higher mass gain from year to year in this short period with b _a = +0.017 m w.e. a⁻², because the increase in the winter balance comes over the summer budget.

3. 3. After the summer balance shift down in τ _s = 1982, the glacier loses mass, with a mean balance of −1.36 m w.e. a⁻¹. Furthermore, due to the concordance of unfavourable trends in winter and summer balances, there is a well-marked trend towards higher mass loss from year to year, with b _a = −0.022 m w.e. a⁻².

These results are in good agreement with the available mean annual geodetic balances mentioned in section 2. Indeed, these independent data indicate that the mean annual balances over the periods 1952–81 and 1981–2003 are −0.36 and −0.99 m a⁻¹, respectively, thus confirming the large increase in mass loss between the two periods.

5.1. Model efficiency and usefulness of its different components

The change-point model M₂ previously discussed introduces two different ruptures for winter and summer balances, which have been shown to be acceptable for the Sarennes series. We have tested a simpler model, M₃, with a single change point, τ, in the two series. The main mode in the distribution of τ is then in 1982 (Fig. 7c), i.e. the year obtained for τ _s with model M₂ (Fig. 7a). This year is identified since it is the most significant of the (τ _w, τ _s) pair. The change-point selection is, however, less clear than with model M₂ because τ = 1982 does not correspond at all to the winter balance series. As a result, the winter trend, w ^*(t), is then very close to zero whatever t (not shown). A single change point must therefore be rejected on the basis of Sarennes data.

We have also analysed the effect of setting a change point τ _a in the annual balance instead of the change point τ _s in the summer balance. This results in a well-marked main mode in 1984 in τ _a’s posterior PDF, thus slightly later than τ _s (Fig. 7d). However, the main difference between this model M₄ and model M₂ is in the trends s ^*(t) and a ^*(t) around the change points (not shown). The mass-balance trend is positive (b _a1 = +0.027 m w.e.a⁻²) until 1984, then becomes slightly negative (b _a2 = −0.013 m w.e. a⁻²). On the other hand, for summer balance, since τ _w ≠ τ _a, three subperiods appear. Having two change points in summer balance and only one in annual balance is highly unlikely and justifies preferring model M₂.

Finally, nonlinear effects have been neglected until here. Removing this assumption requires adding, with suitable constraints to ensure model identifiability, interaction terms , for the annual balance and for the winter balance in model M₂, leading to the full bivariate ANOVA model M₅:

(33)

Among all , only is significantly nonzero at the 95% credibility level for Sarennes. This indicates that only stake 1 behaves differently from the glacier’s mean in some specific years, and only for winter balance. This can be explained by enhanced winter accumulation due to snowdrifts/ avalanches in the area of stake 1. As discussed in section 4, these effects were part of the residuals with models M₁–M₄ and are, with this model, captured in the interaction term.

Second-order temporal terms (δ_t ’s) are all very small, indicating no year with excess or deficits in mass balance and/or winter balance that are not captured by the β_t ’s (Fig. 9). This confirms that most of the temporal structure is actually captured by the linear temporal terms. Note that the δ_t ’s are somewhat more significant for winter than for annual balance, which makes sense since it has been established with model M₁ that the decomposition between space and time is more adapted to the data for annual balance than for winter balance.

Fig. 9. Box plots of the additional nonlinear temporal terms in the annual (a) and winter (b) balances, model M₅. The posterior mean is shown in black, the interquartile range in grey and the 95% credibility interval by dashed lines. Temporal cross terms and are displayed in physical units, whereas spatial terms, γ, are arbitrarily normalized to unity.

Adequacy between model and data measured in terms of variance of residuals or using a model selection tool such as the deviance information criterion (DIC; Reference Spiegelhalter, Best, Carlin and van der LindeSpiegelhalter and others, 2002) is better with model M₅ than with model M₂ (and than with any other model) because model M₅ is much more flexible due to the new parameters introduced. However, standard deviations for the three interaction terms (estimated, e.g., as for annual balance) range from 0.122 to 0.185 m w.e. a⁻¹only for the three components of the balance. This is not negligible, but not high enough to represent a large part of data variability, so that the normality of the residuals is still rejected with model M₅ for winter balance (Kolmogorov–Smirnov test). Furthermore, as a consequence of model structure, the three obtained β_t series and underlying trends are fully identical with model M₅ to those provided by model M₂. Model M₅ cannot therefore be considered as really superior to model M₂, being much more complex in terms of specification and inference without providing much valuable additional information.

5.2. Interpretation of temporal patterns

In this work, fluctuations of different seasonal components in mass balance have been related to time taken as an explicit covariate rather than to seasonal climatic descriptors such as cumulated winter snow depth or summer mean temperature. The main justification of such a time-explicit framework is its flexibility: it does not impose a forced link with climate, and therefore allows the extraction of refined temporal patterns as close as those present in the data. On the other hand, their climatic relevance is then not granted, and one must be extremely cautious when considering trends, change points, and their interpretation in term of climatic drivers.

For instance, Sarennes has a cumulative wide balance of −48 m w.e. over the period of record, and there is an average elevation decrease of 53 m over the entire glacier surface. We would therefore expect an enhanced surface energy exchange due to the elevation gradient of mass balance. Assuming a value of −0.66 m w.e. a⁻¹(100 m)⁻¹ for the mass-balance gradient, a mean annual elevation decrease of −0.89 m a⁻¹ over six decades results in a negative trend in mass balance of −0.0059 m w.e. a⁻² due to the lowering of the glacier surface. This represents nearly 27% of the trend detected by model M₂ in the annual balance time series, which is not a direct consequence of changes in climatic forcing but rather a feedback from the glacier. Another potential source of glacier topography retroaction is the decrease in area at low elevations which removes negative contributions to the glacier-total balance (Reference Elsberg, Harrison, Echelmeyer and KrimmelElsberg and others, 2001). However, this effect has no impact here because the terms are not spatial integrated balances, but uniform annual deviations from the mean surface balance at all individual measurement locations.

Nevertheless, other observations support the assumption that the major explanation of our results is recent climate change. For example, the same trends in summer balances have been reported for other French and Swiss glaciers with analogous temporal break points (Reference VincentVincent, 2002; Reference Huss, Bauder, Funk and HockHuss and others, 2008). The rise in winter precipitation since 1977 that has been found has a more regional signification. Indeed, no specific trend is reported in glacier winter balances in the eastern Swiss Alps, for instance on Silvretta and Clariden glaciers (Reference Müller and KappenbergerMüller and Kappenberger, 1991; Reference Huss and BauderHuss and Bauder, 2009). Nevertheless, closer to Sarennes, Aletschgletscher, Valais, Switzerland, displays a synchronous change point in winter balance around 1974–76, though less significant in magnitude than for Sarennes (Reference Huss and BauderHuss and Bauder, 2009). Furthermore, this change point is present in various snow and weather series in the French Alps (Reference Durand, Laternser, Giraud, Etchevers, Lesaffre and MérindolDurand and others, 2009a,Reference Durand, Giraud, Laternser, Etchevers, Mérindol and Lesaffreb) and avalanche runout elevation series (Reference Eckert, Parent, Kies and BayaEckert and others, 2010c), which supports its climatic relevance, at least at the scale of the French Alps and their close surroundings.

These preliminary interpretations have now to be expanded by analysing in detail the correlations between the different mass-balance series obtained, their underlying trends, and seasonal climatic data. This was not done in this study because its scope was mainly methodological. The different models’ components have therefore been used much more to discuss the advantages of our approach compared to previous ones than to evaluate the contribution of Glacier de Sarennes as a proxy indicator of climate change around the French Alps. Hence, this will be the scope of further work which will involve cross-comparison of the refined indicators provided by our analysis of Sarennes’ exceptionally long and complete series with systematic variations of the constraining climatic factors at large and local scales.

5.3. Implication for the overall and cumulative balances

In a previous paper (Reference Thibert and VincentThibert and Vincent, 2009), a combination of the variance decomposition of the annual mass-balance series using Lliboutry’s classical linear model with aerial photogrammetry was proposed. By minimizing systematic errors, this gives the best possible estimate of the glacier-wide balance. The same goal can be achieved with model M₂ instead of Lliboutry’s model over the periods 1952–81 and 1981–2003 on which, as mentioned in section 2, the mean geodetic balances are available. Note that even if the period of record is only divided into two periods, before and after 1981, instead of the three periods identified by model M₂ (before τ _w = 1976, between τ _w and τ _s = 1982, and after τ _f), 1981 is approximately the middle of the 6–7year long 1977–84 period of high mass-balance regime.

As indicated by Reference Thibert and VincentThibert and Vincent (2009), because of the linearity of the variance decomposition, the glacier-wide balance, , can be evaluated from the annual deviation, , and the mean annual balance,Φ 4), obtained from photogrammetry according to Figure 2b:

(34)

The upper-case letters designate spatially integrated values as opposed to point measurements such as in Equation (2), for example. In our case, Φ takes two different values, before and after 1981. Furthermore, since no photogrammetry measurement is available for the initial and final years of record, the glacier total balance is calculated from the terms. The cumulative balance, , for glaciological year t is then given by:

(35)

Under model M₂, is centred on a linear time trend with slope b _a. which is homogeneous to a rate. The sum of the terms then has a second- order time dependence such as:

(36)

where C is a constant which can be determined for each of the two considered periods using the control mass balances from photogrammetry and by imposing continuity around the change points. Over each period, the cumulative glacier-wide balance is therefore a parabolic function of time, with curvature given by the sum in winter and summer balance trends:

(37)

Considering the two periods given by photogrammetry instead of the three given by the variance analysis is thus equivalent to neglecting the slight increase in the mass-balance trend between 1977 and 1982 (from +0.011 to +0.017 m w.e. a⁻²).

The obtained cumulative glacier-wide balance perfectly matches the parabolic trend in 1952, 1981 and 2003 as a result of the fitting of the mean annual balance to photogrammetry measurements performed on these dates (Fig. 10). For all other years, the obtained parabolic curve also fits the data well, which is remarkable since it results from integrating the mean trend provided by model M₂ rather than from direct fitting to cumulated data at each stake. This provides another argument in favour of the assumption of robust underlying linear trends for the changes in mean winter and summer balances.

Fig. 10. Cumulative glacier-wide balance, , obtained by combining photogrammetry performed in 1952, 1981 and 2003, and annual terms centred on a linear time trend given by model M₂. As a result of the two linear trends in annual balances, trends in cumulative balance, , are parabolic over the 1949–81 and 1981–2007 periods.

For instance, the curvature of the cumulative balance is positive before 1981 (+0.011 m w.e. a⁻²) and negative since (−0.022 m w.e. a⁻²). In mean, i.e. when the predominant temporal pattern is separated from interannual variability modelled as random, mass loss decelerated before 1981 while Sarennes was approaching an equilibrium state, which was reached in the early 1980s. On the other hand, in the last 25 years, the predominant pattern is an accelerated glacier mass loss clearly indicated by the negative curvature of the cumulative balance. This shows that Sarennes is rapidly moving away from the last equilibrium state of the 1980s. Even considering that a part of this acceleration in mass loss (27%) is retroactive due to glacier topography, this confirms that the glacier is now still far from its equilibrium state under the current climate and/ or that climate change is still in progress.