# Extracting the temporal signal from a winter and summer mass-balance series: application to a six-decade record at Glacier de Sarennes, French Alps

Nicolas Eckert, H. Baya, Emmanuel Thibert, C. Vincent

## ► To cite this version:

Nicolas Eckert, H. Baya, Emmanuel Thibert, C. Vincent. Extracting the temporal signal from a winter and summer mass-balance series: application to a six-decade record at Glacier de Sarennes, French Alps. Journal of Glaciology, 2011, 57 (201), pp.134-150. 10.3189/002214311795306673 . insu-00604753

HAL Id: insu-00604753

https://insu.hal.science/insu-00604753

Submitted on 25 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Extracting the temporal signal from a winter and summer mass-balance series: application to a six-decade record at Glacier de Sarennes, French Alps

N. ECKERT,<fnref n="1" /> H. BAYA,<fnref n="1" /> E. THIBERT,<fnref n="1" /> C. VINCENT<fnref n="2" />

**ABSTRACT.** Temporal trends related to recent climatic fluctuations are extracted from the longest glacier-wide winter and summer mass-balance series recorded in the Alps, at Glacier de Sarennes, France. For this, all point balances measured at the glacier surface are used, and different statistical models are developed and tested. First, Lliboutry's linear variance analysis model is extended to the two seasonal components of the balance. The explicit modelling of variability sources and correlations is proved useful for appropriately quantifying uncertainties in the different components of the balance and estimating missing data. Next, a non-exchangeable structure is added to model the winter and summer balance time series. Two change points separating different underlying trends are thus detected. The first change was in 1976, with a shift of +23% in the winter balance. The second was in 1982 for the summer balance series. These systematic changes explain 20-30% of the variability of the different components of the balance, the rest being made up of random interannual fluctuations. Simplified and/or less physically based models are less efficient in capturing data variability. As a result, the cumulative glacier-wide balance shows systematic parabolic trends, which result in an accelerated mass loss for Glacier de Sarennes over the last 25 years.

## 1. INTRODUCTION

In the European Alps, glacier mass balance results mainly from accumulation (mass input) of snow in winter and ablation (mass loss) due to the melting of snow and ice in summer. Fluctuations of mass-balance seasonal components are therefore strongly related to climate fluctuations. More generally, glacier mass balances are well-known useful indicators of climate change over the last few centuries (Oerlemans and Fortuin, 1992; Haeberli, 1995; Vincent, 2002; Vincent and others, 2004; Kaser and others, 2006; Solomon and others, 2007).

However, inferring the climatic signal from a glacier mass-balance series is not an easy task. For example, the glacier-wide balance obtained from spatial averaging of point mass balances (measured at individual locations) is a function of glacier geometry (Ohmura and others, 2007). To overcome this problem, Elsberg and others (2001) have proposed a method that corrects surface and elevation changes in the integrated glacier-wide balance. An alternative is to directly use surface point mass balances (Rasmussen, 2004; Vincent and others, 2004; Huss and Bauder, 2009; Thibert and Vincent, 2009), but this implies separating the overall annual effect from the spatial variability at the glacier surface. This can be done using variance decompositions (analysis of variance: ANOVA).

The suitability of annual point mass-balance time series for linear variance decomposition was first proposed by Meier and Tangborn (1965) and Hoinkes (1970). Lliboutry (1974) generalized the approach by proposing what is formally an incomplete two-way ANOVA model:

$$b_{it}^{a} \sim N_{1}(\alpha_{i}^{a} + \beta_{t}^{a}, \sigma_{a}^{2}), \quad (1)$$

where $N_1$ denotes the normal distribution and $b_{it}^a$ the annual

(superscript a) point mass balance observed at spatial coordinates $(x_i, y_i)$ during the glaciological year, $t \in [t_0, t_0 + T - 1]$, where $t_0$ is the first year considered and $T$ is the total number of years with observations. $\alpha_i^a = \alpha^a(x_i, y_i)$ is the spatial term. The notation indicates that the measurements are made using stakes $i \in [1, N]$, for which the number and positions are chosen to provide a good vision of the spatial variability at the glacier surface. The constraint $\sum_{t=t_0}^{t_0+T-1} \beta_t^a = 0$ is used to ensure model identifiability, so that $\beta_t^a$ is the annual deviation. The residuals of standard deviation, $\sigma_a^2$, are assumed to be independent and Gaussian. The model is referred to as incomplete because the interaction terms between space and time are not distinguished from the residuals.

Evaluating the glacier-wide balance involves integration over the glacier surface. Under an assumption of equally weighted measurement points, the glacier-wide mean balance over the studied period is $\langle B^a \rangle = \frac{1}{N} \sum_{i=1}^{N} \alpha_i^a$. Note that the way spatial averaging is performed has little influence on the annual deviations, and that implementation of other choices (e.g. area-altitude distribution weighting coefficients) than the direct arithmetic mean used here is straightforward. Furthermore, the annual glacier-wide balance is then simply $B_t^a = \langle B^a \rangle + \beta_t^a$, upper-case letters denoting glacier-wide variables. This highlights that the mean annual deviation of the balance $\beta_t^a$ at all measurement sites equals the glacier-wide annual deviation. Therefore, the $\beta_t^a$ time series extracted by the model from point balances also apply to the wide balance, and can be used to characterize the response of the glacier to climate fluctuations.

Because of its simplicity and ability to fit various point mass-balance series, this model has met with great success in

<footnote n="1">¹Cemagref, UR ETGR, 2 rue de la Papeterie, BP 76, 38402 Saint-Martin-d'Hères Cedex, France
E-mail: emmanuel.thibert@cemagref.fr</footnote>

<footnote n="2">²Laboratoire de Glaciologie et Géophysique de l'Environnement, CNRS/Université Joseph Fourier - Grenoble I,
54 rue Molière, BP 96, 38402 Saint-Martin-d'Hères Cedex, France</footnote>
the glaciological community (Vallon and Leiva, 1981; Kuhn, 1984; Rasmussen, 2004; Thibert and others, 2008; Thibert and Vincent, 2009). However, some studies have reported discrepancies between model results and real datasets, caused by unsteady spatial terms or temporal terms decreasing at high elevations (Oerlemans and Hoogendoorn, 1989; Funk and others 1997; Vallon and others, 1998; Dyurgerov and Dwyer, 2001; Soruco and others, 2009).

From a more methodological point of view, a first limitation of the model proposed by Lliboutry (1974) is that it does not explicitly consider the two specific contributions of summer and winter balance components to the annual balance. Rasmussen and Andreassen (2005) therefore applied their variance analysis separately to both winter and summer components of the balance, with good agreement between model results and data for studies concerning Norway. Although it represents a big step forward from a physical point of view, such an approach remains questionable. Indeed, by modelling the different components of the balance separately, it neglects the fact that they are correlated. Their correlations must therefore be estimated empirically from the obtained annual terms. In statistics, empirical estimation denotes a procedure that does not explicitly rely on a model fitted on data using a procedure such as mean square error minimization. This leads to different estimates, with generally higher estimation errors. Furthermore, among winter, summer and annual balance, only two contributions are usually measured (the third being calculated from), so that neglecting the correlation between the two measured variables may imply underestimating the uncertainty around the third unobserved variable.

A second limitation of Lliboutry's model is that it implicitly considers stationary conditions over the period of record, with random annual fluctuations around a mean balance, which is inappropriate in the context of a marked climate change. Standard time series can be used to extract structured patterns (e.g. memory effects, systematic trends, break points) from the $β_t$ series, with the ultimate goal of relating them to similar behaviours in more direct climatic data series. However, this poses an even greater statistical problem. Indeed the $β_t$'s are unobserved, so that the uncertainty related to their estimation must be taken into account when estimating the investigated temporal patterns. Not doing so may lead to an overestimation of the significance levels of the temporal patterns because the studied series are only a few decades long.

The objective of this paper is therefore to reconsider the problem of extracting a climatic trend from a glacier mass balance recorded by a sampling site network using a variance decomposition. To do so, we first propose an expansion of Lliboutry's model adapted to correlated pairs. Second, we show how hierarchy (Wikle, 2003; Clark and Gelfand, 2006) can be used to isolate the annual terms and explicitly model their temporal structure within the same statistical framework (Carlin and others, 1992; Banerjee and others, 2004), so as to relax in a consistent manner the underlying assumption of stationary conditions. Inspired by work in other fields (Clark, 2005; Eckert and others, 2007, 2010c), implementation is performed using Bayesian schemes thanks to the availability of simulation-based methods for parameter estimation (Tanner, 1992; Gaucherel and others, 2008). The proposed approach and for instance the effectiveness of the different model components are illustrated using the six-decade Sarenes mass-balance series.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Location of stakes where balances are measured on Glacier de Sarenes. Contours on the glacier are at 25 m intervals.</figcaption></figure>

## 2. DATA AND SITE

The mass-balance series of Glacier de Sarenes (45°07′N, 6°07′E), Massif des Grandes Rousses, France, has been measured since $t_0=1949$ (Valla, 1989; Thibert and others, 2008). Sarenes is a small south-facing glacier whose surface decreased from $0.85\ \text{km}^2$ in 1952 to $0.41\ \text{km}^2$ in 2003. Over the period of record, the elevation ranges from 2800 to 3100 m a.s.l. (Fig. 1). On this glacier, the winter balance, $b_{it}^w$, is measured from cores drilled at the end of the winter. The annual balance, $b_{it}^a$, is measured in late summer from emergence variations of stakes inserted in the ice of the ablation zone and cores drilled in the accumulation zone. The summer balance term at each location, $b_{it}^s$, is then estimated from $b_{it}^s = b_{it}^a - b_{it}^w$. Point balances ($b_{it}^a$, $b_{it}^w$) are measured at five sites for which the mean elevation over the period of record is respectively 2860, 2920, 2945, 2995 and 3020 m a.s.l. Mean altitude change over the period of record is –36 m, ranging from –61 m in the lower part, to –15 m in high-elevation areas. From 1949 to 1957, only two to four sites were sampled, and only four since 2003 because of the retreat of the tongue. The 1949–2008 experimental table is therefore a partially complete dataset for the three terms of the balance (27 missing values). These annual data are expressed in m w.e. $\text{a}^{-1}$ and are published on http://www.lgge.ujf-grenoble.fr/ServiceObs/.

Measurements are not performed at fixed dates: the investigated annual mass changes refer to the so-called floating date system. Durations of mass-balance year are therefore variable, and account for both the effects of changes in accumulation and ablation rates, and of their duration. Note that, as field visits do not coincide in general with the glacier minimum or maximum mass, data differ somewhat from the theoretical stratigraphic system in which annual balance is defined between one year's minimum and the next year's, and winter balance between one year's maximum and the previous year's. However, we assume that these differences are small, and do not apply a correction.

A routinely calculated glacier-wide balance (published in World Glacier Monitoring Service (WGMS) glacier mass-balance bulletins) uses an integration relationship based on area-altitude distribution weighting coefficients (Valla, 1989; Thibert and others, 2008). Over the 1949–2007 period, this calculation gives a mean annual glacier-wide
balance of -0.92 m w.e. a⁻¹, while mean annual glacier-wide winter and summer balances are 1.68 and -2.60 m w.e. a⁻¹, respectively. Additional annual mass-balance data estimated by photogrammetry (geodetic balance) are available in 1981 (Valla and Piedallu, 1997), 1952 and 2003 (Thibert and others, 2008) for comparison and control.

# 3. METHODS

## 3.1. Bivariate formulation of Lliboutry's model

To generalize Lliboutry's model to the actually observed pair $(b_{it}^a, b_{it}^w)$, we introduce spatial terms $\alpha_i^w$ and temporal terms $\beta_t^w$ for winter balance under the additional constraint $\sum_{t=t_0}^{t_0+T-1} \beta_t^w = 0$:

$$ \begin{pmatrix} b_{it}^a \\ b_{it}^w \end{pmatrix} \sim N_2 \left[ \begin{pmatrix} \alpha_i^a + \beta_t^a \\ \alpha_i^w + \beta_t^w \end{pmatrix}, \begin{pmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{pmatrix} \right], \quad (2) $$

where superscript w denotes the winter period. The overall mean $\langle B^w \rangle = \frac{1}{N} \sum_{i=1}^N \alpha_i^w$ is the mean glacier-wide winter balance over the studied period, whereas the $\beta_t^w$ series represents the annual deviations to the mean. $N_2$ denotes the two-dimensional (2-D) Gaussian vector of variance-covariance matrix $\Sigma$ which includes the marginal variances $\Sigma_{11} = \text{Var}(b_{it}^a) = \sigma_a^2$ and $\Sigma_{22} = \text{Var}(b_{it}^w) = \sigma_w^2$. The residual pair $(\varepsilon_{it}^a, \varepsilon_{it}^w)$ 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:

$$ \left\{ \begin{aligned} \varepsilon_{it}^a &= b_{it}^a - (\alpha_i^a + \beta_t^a) \\ \varepsilon_{it}^w &= b_{it}^w - (\alpha_i^w + \beta_t^w) \end{aligned} \right. . \quad (3) $$

The main advantage of the joint formulation is that the covariance $\Sigma_{12} = \Sigma_{21} = \text{Cov}(b_{it}^a, b_{it}^w)$ is also available, providing easy access to the correlation coefficient $\rho_{b^a b^w}$ between the residual pair $(\varepsilon_{it}^a, \varepsilon_{it}^w)$:

$$ \rho_{b^a b^w} = \frac{\text{Cov}(b_{it}^a, b_{it}^w)}{\sqrt{\text{Var}(b_{it}^a)} \sqrt{\text{Var}(b_{it}^w)}} = \frac{\Sigma_{12}}{\sigma_a \sigma_w}. \quad (4) $$

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

The pair of annual deviations $(\beta_t^a, \beta_t^w)$ is assumed to be sampled from a 2-D Gaussian vector of variance-covariance matrix $\Sigma_\beta$. Its mean is set to zero for consistency with the imposed constraints:

$$ \left( \begin{pmatrix} \beta_t^a \\ \beta_t^w \end{pmatrix} \right) \sim N_2 \left[ \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \Sigma_\beta \right]. \quad (5) $$

The hidden quantities $(\beta_t^a, \beta_t^w)$ 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 $\Sigma_\beta$ covariance matrix. The marginal variances $\Sigma_{\beta 11} = \text{Var}(\beta_t^a) = \sigma_a^2$ and $\Sigma_{\beta 22} = \text{Var}(\beta_t^w) = \sigma_w^2$ 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 $\Sigma_{\beta 21} = \Sigma_{\beta 12} = \text{Cov}(\beta_t^a, \beta_t^w)$, which provides easy access to their correlation coefficient $\rho_{\beta^a \beta^w}$. It is definitely non-negligible, since mass balance physically

depends on winter balance:

$$ \rho_{\beta^a \beta^w} = \frac{\text{Cov}(\beta_t^a, \beta_t^w)}{\sqrt{\text{Var}(\beta_t^a)} \sqrt{\text{Var}(\beta_t^w)}} = \frac{\Sigma_{\beta 12}}{\sigma_{\beta^a} \sigma_{\beta^w}}. \quad (6) $$

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

$$ b_{it}^s = b_{it}^a - b_{it}^w, \quad (7) $$

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

$$ b_{it}^s \sim N_1 (\alpha_i^s + \beta_t^s, \sigma_s^2), \quad (8) $$

and its variance decomposition into a spatial term $\alpha_i^s$, an annual term $\beta_t^s$ verifying $\sum_{t=t_0}^{t_0+T-1} \beta_t^s = 0$ and a residual $\varepsilon_{it}^s$ is readily available:

$$ \left\{ \begin{aligned} \alpha_i^s &= \alpha_i^a - \alpha_i^w \\ \beta_t^s &= \beta_t^a - \beta_t^w \\ \varepsilon_{it}^s &= b_{it}^s - (\alpha_i^s - \beta_t^s) = \varepsilon_{it}^a - \varepsilon_{it}^w \end{aligned} \right. . \quad (9) $$

Keeping an arithmetic averaging over the glacier surface, the mean summer balance over the studied period is then simply $\langle B^s \rangle = \frac{1}{N} \sum_{i=1}^N \alpha_i^s = \langle B^a \rangle - \langle B^w \rangle$. Moreover, marginal variances for the summer balance residual $\varepsilon_{it}^s$ and for the annual deviation to the mean summer balance $\beta_t^s$ can be computed from the covariance matrices $\Sigma$ and $\Sigma_\beta$, 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 $\sigma_a + \sigma_w$ for $\sigma_s$, and smaller than $\sigma_{\beta^a} + \sigma_{\beta^w}$ for $\sigma_{\beta^s}$:

$$ \sigma_s^2 = \Sigma_{11} + \Sigma_{22} + 2\Sigma_{12}, \quad (10) $$

$$ \sigma_{\beta^s}^2 = \Sigma_{\beta 11} + \Sigma_{\beta 22} + 2\Sigma_{\beta 12}, \quad (11) $$

$$ \sqrt{\sigma_a^2 + \sigma_w^2} \le \sigma_s \le \sigma_a + \sigma_w, \quad (12) $$

$$ \sqrt{\sigma_{\beta^a}^2 + \sigma_{\beta^w}^2} \le \sigma_{\beta^s} \le \sigma_{\beta^a} + \sigma_{\beta^w}. \quad (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 $\rho_{b^a b^s}$ and $\rho_{b^w b^s}$ between the residuals pairs $(\varepsilon_{it}^a, \varepsilon_{it}^s)$ and $(\varepsilon_{it}^w, \varepsilon_{it}^s)$ are given by:

$$ \rho_{b^a b^s} = \frac{\text{Cov}(b_{it}^a, b_{it}^s)}{\sqrt{\text{Var}(b_{it}^a)} \sqrt{\text{Var}(b_{it}^s)}} = \frac{1}{N^T} \sum_{i=1}^N \sum_{t=t_0}^{t_0+T-1} \frac{\varepsilon_{it}^a \varepsilon_{it}^s}{\sigma_a \sigma_s}, \quad (14) $$

$$ \rho_{b^w b^s} = \frac{1}{N^T} \sum_{i=1}^{N} \sum_{t=t_0}^{t_0+T-1} \frac{\varepsilon_{it}^w \varepsilon_{it}^s}{\sigma_w \sigma_s}, \quad (15) $$

whereas the correlation coefficients $\rho_{\beta^a \beta^s}$ and $\rho_{\beta^w \beta^s}$ between the annual deviations are given by:

$$ \rho_{\beta^a \beta^s} = \frac{\text{Cov}(\beta_t^a, \beta_t^s)}{\sqrt{\text{Var}(\beta_t^a)} \sqrt{\text{Var}(\beta_t^s)}} = \frac{1}{T} \sum_{t=t_0}^{t_0+T-1} \frac{\beta_t^a \beta_t^s}{\sigma_{\beta^a} \sigma_{\beta^s}}, \quad (16) $$

$$ \rho_{\beta^w \beta^s} = \frac{1}{T} \sum_{t=t_0}^{t_0+T-1} \frac{\beta_t^w \beta_t^s}{\sigma_{\beta^w} \sigma_{\beta^s}}. \quad (17) $$

Equation (2) separates different independent effects using variance decomposition. We introduce the ratios $\lambda_a, \lambda_w$ and $\lambda_s$ (Equations (18–20)) to quantify the specific contribution
<figure><img src="image_3.png" /><figcaption>**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 αᵢ's spatial mean ⟨Bᵃ⟩ is also indicated.</figcaption></figure>

(% 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 αᵢ terms. This is reasonable for a small glacier such as Sarenes (Thibert and Vincent, 2009), but would not be appropriate for a wider glacier covering a large elevation range:

$$ \lambda_a = \frac{\sigma_{\beta a}^2}{\sigma_{\beta a}^2 + \sigma_a^2}, \qquad (18) $$

$$ \lambda_w = \frac{\sigma_{\beta w}^2}{\sigma_{\beta w}^2 + \sigma_w^2}, \qquad (19) $$

$$ \lambda_s = \frac{\sigma_{\beta s}^2}{\sigma_{\beta s}^2 + \sigma_s^2}. \qquad (20) $$

Figure 2a summarizes this first model, noted M₁, using a direct acyclic graph (Lauritzen, 1996). On the one hand, the parameters (Σ, Σ_β, α_i^a, α_i^w), i ∈ [1, N] and hidden variables (β_t^a, β_t^w), t ∈ [t₀, t₀ + T - 1] 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 (Σ, Σ_β, α_i^a, α_i^w) and (β_t^a, β_t^w), which somewhat restricts the dimension of the estimation problem. With regard to Rasmussen 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.

## 3.2. Extracting structured temporal trends

With model M₁, the annual terms (β_t^a, β_t^w) 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 βᵢ's. Moreover, for better consistency with the physical processes involved, we prefer working with the seasonal pair (β_t^w, β_t^s) rather than with the pair (β_t^a, β_t^w).

A time-series analysis of the (β_t^w, β_t^s)'s can theoretically be achieved using various time-series models developed in other fields, such as econometrics (Booth and Smith, 1982), hydrology (Rao and Titrojondro, 1996), meteorology (Bloomfield, 1992) or climatology (von 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 (Mearns and others, 1997), variance (Diaz, 1982) or extreme values (Nogaj 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 Sarenes data (see section 4) and has already been applied to similar problems in hydrology (Perreault and others, 2000a,b).

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 β_t^w's (Equations (21) and (22)) and the β_t^s's (Equations (23) and (24)) are broken down into random noises and structured linear trends. The random noises, with respective variances σ_w^2 and σ_s^2, model the unstructured interannual fluctuations. The notation σ_s^2, for example, indicates that the parameter σ_s^2 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:

$$ \beta_t^w \sim N_1(a_{w1} + b_{w1} t, \sigma_{\beta w_1}^2), t \in [t_0, \tau_w], \quad (21) $$

$$ \beta_t^w \sim N_1(a_{w2} + b_{w2} t, \sigma_{\beta w_2}^2), t \in [\tau_w + 1, t_0 + T - 1], \quad (22) $$

$$ \beta_t^s \sim N_1(a_{s1} + b_{s1} t, \sigma_{\beta s_1}^2), t \in [t_0, \tau_s], \quad (23) $$

$$ \beta_t^s \sim N_1(a_{s2} + b_{s2} t, \sigma_{\beta s_2}^2), t \in [\tau_s + 1, t_0 + T - 1]. \quad (24) $$

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

$$ a^*(t) = w^*(t) + s^*(t) = (a_w + a_s) + (b_w + b_s)t. \quad (25) $$

If $\tau_s = \tau_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 $\tau_s \neq \tau_w$, three possibly different linear trends and intercepts are obtained for respectively $t \le \min(\tau_s, \tau_w)$, $t \in [\min(\tau_s, \tau_w) + 1, \max(\tau_s, \tau_w)]$ and $t \ge \max(\tau_s, \tau_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₁:

$$ \sigma_{\beta^a}^2 = E(\beta_t^{a2}) = \frac{1}{T} \sum_{t=t_0}^{t_0+T-1} (\beta_t^a)^2, \quad (26) $$

$$ \sigma_{\beta^w}^2 = E(\beta_t^{w2}) = \frac{1}{T} \sum_{t=t_0}^{t_0+T-1} (\beta_t^w)^2, \quad (27) $$

$$ \sigma_{\beta^s}^2 = E(\beta_t^{s2}) = \frac{1}{T} \sum_{t=t_0}^{t_0+T-1} (\beta_t^s)^2. \quad (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 co-efficient correlations $\rho_{\beta^a \beta^a}$ and $\rho_{\beta^w \beta^w}$ can be evaluated using Equations (16) and (17) respectively. For $\rho_{\beta^a \beta^w}$, 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₂:

$$ \rho_{\beta^a \beta^w} = \frac{\frac{1}{T} \sum_{t=t_0}^{t_0+T-1} \beta_t^a \beta_t^w}{\sigma_{\beta^a} \sigma_{\beta^w}}. \quad (29) $$

Finally, to compare the respective contributions to the annual variability of the temporal structure and the random fluctuations, we define the global ratios $\mu_{ar}$, $\mu_w$ and $\mu_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^2$ in a

non-hierarchical linear regression:

$$ \mu_a = \frac{\mathrm{Var}(a^*(t))}{\mathrm{Var}(a^*(t)) + \sigma_{\beta^a}^2}, \quad (30) $$

$$ \mu_w = \frac{\mathrm{Var}(w^*(t))}{\mathrm{Var}(w^*(t)) + \sigma_{\beta^w}^2}, \quad (31) $$

$$ \mu_s = \frac{\mathrm{Var}(s^*(t))}{\mathrm{Var}(s^*(t)) + \sigma_{\beta^s}^2}. \quad (32) $$

Obviously, the different proposed ratios as well as the correlation coefficients between the $\beta_t$'s can also be computed separately before and after the change points $\tau_w$ and $\tau_s$ instead of globally over the entire observation period. This is, for instance, straightforward for $\lambda_{s}$, since it only implies computing the ratio with the two values $\sigma_{\beta^1}$ and $\sigma_{\beta^2}$. 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 ($a_{w,}$, $a_{s,}$, $b_{w,}$, $b_{s,}$, $\sigma_{\beta^w,}$, $\sigma_{\beta^a,}$, $\tau_w$, $\tau_s$) 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 $\Sigma_{\beta12}$.

### 3.3. Bayesian inference and associated MCMC sampling

The different models proposed in this paper have been inferred within a Bayesian framework (Berger, 1985; Parent 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 Gilks and others (1996) and Brooks (1998). Examples of applications in a related field are given by Eckert and others (2008, 2010b). In our case, the efficient Gibbs sampler (Geman 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 (Brooks 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
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

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 $\sigma_a$ and $\sigma_{\beta^a}$ are direct results of the MCMC procedure. The others were obtained by computing their values at each step of the sampling procedure.

# 4. RESULTS

## 4.1. Lliboutry's bivariate model

Table 1 summarizes the estimates obtained by applying model M₁ to the Sareennes series. It appears that $\langle B^a \rangle$, 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 $\sim 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 ($\langle B_{AAD}^a \rangle$, $\langle B_{AAD}^w \rangle$, $\langle B_{AAD}^S \rangle$; Table 1). They are also clearly consistent with previous studies of the same glacier (Valla, 1989; Vincent and Vallon, 1997; Thibert and

others, 2008; Thibert and Vincent, 2009) and other glaciers in the French and Swiss Alps (Müller and Kappenberger, 1991; Vincent, 2002; Vincent and others, 2004; Huss and Bauder, 2009).

Marginal variances ($\sigma_{a_r}$, $\sigma_{w_r}$, $\sigma_{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 $\sigma_{\beta^a}$, $\sigma_{\beta^w}$, $\sigma_{\beta^S}$ 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 $\alpha_i$ series, i.e. not more than 0.044 m w.e. a⁻¹. All this indicates that the interannual variability between the $\beta_i$'s is the major source of variability at the glacier scale, constraining the three ratios $\lambda_{a_r}$, $\lambda_{w_r}$ and $\lambda_{s}$ to be above 90%. It also a posteriori justifies neglecting the spatial variability while computing the ratios of Equations (18–20).

Residuals account for nonlinear effects that are not systematic in space and/or time and cannot therefore be captured in the $\alpha$ and $\beta$ 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 $\rho_{b^a b^w}$ and $\rho_{b^a b^S}$ are easily understandable since when a local excess of winter mass
<table><caption>**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</caption><thead><tr><th></th><th>Mean</th><th>SD</th><th>2.5%</th><th>97.5%</th></tr></thead><tbody><tr><td>⟨B<sup>a</sup>⟩ (m w.e. a<sup>-1</sup>)</td><td>-0.97</td><td>0.019</td><td>-1.01</td><td>-0.93</td></tr><tr><td>⟨B<sup>w</sup>⟩ (m w.e. a<sup>-1</sup>)</td><td>1.69</td><td>0.014</td><td>1.67</td><td>1.72</td></tr><tr><td>⟨B<sup>s</sup>⟩ (m w.e. a<sup>-1</sup>)</td><td>-2.66</td><td>0.019</td><td>-2.70</td><td>-2.62</td></tr><tr><td>⟨B<sup>a</sup><sub>AAD</sub>⟩ (m w.e. a<sup>-1</sup>)</td><td>-0.92</td><td></td><td></td><td></td></tr><tr><td>⟨B<sup>w</sup><sub>AAD</sub>⟩ (m w.e. a<sup>-1</sup>)</td><td>1.68</td><td></td><td></td><td></td></tr><tr><td>⟨B<sup>s</sup><sub>AAD</sub>⟩ (m w.e. a<sup>-1</sup>)</td><td>-2.60</td><td></td><td></td><td></td></tr><tr><td>σ<sub>a</sub> (m w.e. a<sup>-1</sup>)</td><td>0.299</td><td>0.015</td><td>0.272</td><td>0.33</td></tr><tr><td>σ<sub>w</sub> (m w.e. a<sup>-1</sup>)</td><td>0.214</td><td>0.011</td><td>0.194</td><td>0.236</td></tr><tr><td>σ<sub>s</sub> (m w.e. a<sup>-1</sup>)</td><td>0.43</td><td>0.021</td><td>0.391</td><td>0.474</td></tr><tr><td>ρ<sub>b<sup>a</sup>b<sup>w</sup></sub></td><td>0.39</td><td>0.06</td><td>0.27</td><td>0.50</td></tr><tr><td>ρ<sub>b<sup>a</sup>b<sup>s</sup></sub></td><td>0.51</td><td>0.05</td><td>0.42</td><td>0.60</td></tr><tr><td>ρ<sub>b<sup>s</sup>b<sup>a</sup></sub></td><td>-0.23</td><td>0.03</td><td>-0.30</td><td>-0.17</td></tr><tr><td>σ<sub>β<sup>a</sup></sub> (m w.e. a<sup>-1</sup>)</td><td>1.123</td><td>0.108</td><td>0.938</td><td>1.357</td></tr><tr><td>σ<sub>β<sup>w</sup></sub> (m w.e. a<sup>-1</sup>)</td><td>0.648</td><td>0.062</td><td>0.539</td><td>0.784</td></tr><tr><td>σ<sub>β<sup>s</sup></sub> (m w.e. a<sup>-1</sup>)</td><td>1.487</td><td>0.143</td><td>1.24</td><td>1.797</td></tr><tr><td>ρ<sub>β<sup>a</sup>β<sup>w</sup></sub></td><td>0.36</td><td>0.11</td><td>0.12</td><td>0.57</td></tr><tr><td>ρ<sub>β<sup>a</sup>β<sup>s</sup></sub></td><td>0.51</td><td>0.09</td><td>0.35</td><td>0.70</td></tr><tr><td>ρ<sub>β<sup>s</sup>β<sup>a</sup></sub></td><td>0.03</td><td>0.02</td><td>0.00</td><td>0.07</td></tr><tr><td>λ<sub>a</sub></td><td>0.93</td><td>0.01</td><td>0.90</td><td>0.96</td></tr><tr><td>λ<sub>w</sub></td><td>0.90</td><td>0.02</td><td>0.86</td><td>0.93</td></tr><tr><td>λ<sub>s</sub></td><td>0.92</td><td>0.02</td><td>0.89</td><td>0.95</td></tr></tbody></table>

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 $\rho_{b^w b^s}$ 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, $\rho_{b^w b^s}$ 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 $\rho_{\beta^a \beta^s}$ exists between annual and summer balance. Identically, annual and winter balances are well correlated ($\rho_{\beta^a \beta^w}$ being significantly positive). Finally, winter and summer balances are uncorrelated with $\rho_{\beta^a \beta^s} \approx 0$. 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 ($\beta_t^a, \beta_t^w, \beta_t^s$) 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. $\langle B^a \rangle$, $\langle B^w \rangle$ and $\langle B^s \rangle$ 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,

<table><caption>**Table 2.** Spatial $\alpha_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⁻¹</caption><thead><tr><th></th><th>Mean</th><th>SD</th><th>2.5%</th><th>97.5%</th></tr></thead><tbody><tr><td>α<sub>1</sub><sup>a</sup></td><td>-1.07</td><td>0.043</td><td>-1.15</td><td>-0.98</td></tr><tr><td>α<sub>2</sub><sup>a</sup></td><td>-1.28</td><td>0.04</td><td>-1.36</td><td>-1.20</td></tr><tr><td>α<sub>3</sub><sup>a</sup></td><td>-1.01</td><td>0.044</td><td>-1.10</td><td>-0.93</td></tr><tr><td>α<sub>4</sub><sup>a</sup></td><td>-0.91</td><td>0.04</td><td>-0.98</td><td>-0.83</td></tr><tr><td>α<sub>5</sub><sup>a</sup></td><td>-0.57</td><td>0.043</td><td>-0.65</td><td>-0.49</td></tr><tr><td>α<sub>1</sub><sup>w</sup></td><td>1.87</td><td>0.031</td><td>1.81</td><td>1.93</td></tr><tr><td>α<sub>2</sub><sup>w</sup></td><td>1.62</td><td>0.029</td><td>1.56</td><td>1.67</td></tr><tr><td>α<sub>3</sub><sup>w</sup></td><td>1.54</td><td>0.031</td><td>1.48</td><td>1.60</td></tr><tr><td>α<sub>4</sub><sup>w</sup></td><td>1.64</td><td>0.028</td><td>1.58</td><td>1.70</td></tr><tr><td>α<sub>5</sub><sup>w</sup></td><td>1.80</td><td>0.031</td><td>1.74</td><td>1.86</td></tr><tr><td>α<sub>1</sub><sup>s</sup></td><td>-2.94</td><td>0.042</td><td>-3.03</td><td>-2.86</td></tr><tr><td>α<sub>2</sub><sup>s</sup></td><td>-2.90</td><td>0.039</td><td>-2.97</td><td>-2.82</td></tr><tr><td>α<sub>3</sub><sup>s</sup></td><td>-2.55</td><td>0.043</td><td>-2.63</td><td>-2.47</td></tr><tr><td>α<sub>4</sub><sup>s</sup></td><td>-2.54</td><td>0.039</td><td>-2.62</td><td>-2.47</td></tr><tr><td>α<sub>5</sub><sup>s</sup></td><td>-2.37</td><td>0.042</td><td>-2.46</td><td>-2.29</td></tr></tbody></table>

glacier-wide balance (i.e. $\beta_t^a + \langle B^a \rangle$), 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.

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 Vincent 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 Vincent 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 $\beta_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
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Box plots of $\beta_t$ terms (annual deviations), model $M_1$. (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_1–M_5$. As shown in Figure 2b, the glacier-wide terms can be obtained by adding to each $\beta_t$ term the interannual mean given in Table 1, i.e. $\langle B^a \rangle$, $\langle B^w \rangle$ and $\langle B^s \rangle$, respectively.</figcaption></figure>

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).

Finally, model adequacy with respect to data can be checked by comparing the residuals ($\epsilon_{it}^a, \epsilon_{it}^w$) with the centred linear model with variance covariance $\Sigma$. 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 $\sigma_a^2$ and $\sigma_w^2$ 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.

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 $\rho_{b^a b^w}$. 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,
<figure><img src="image_6.png" /><figcaption>**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, $b_{3t}^a$, and associated uncertainty for the missing values. (b, c) Posterior distribution of mass balance at stake 3 estimated in (b) 1956 and (c) 1957.</figcaption></figure>

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.

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 ($\beta_t^a, \beta_t^W, \beta_t^S$) explains why the overall variances ($\sigma_{\beta^a}, \sigma_{\beta^W}, \sigma_{\beta^S}$) and correlation coefficients ($\rho_{\beta^a\beta^W}, \rho_{\beta^a\beta^S}, \rho_{\beta^S\beta^W}$) differ slightly with model M₂ (Table 3) from their previous values (Table 1).

More precisely, all marginal variances ($\sigma_{\beta^a}, \sigma_{\beta^W}, \sigma_{\beta^S}$) are reduced, and $\sigma_{\beta^S}$ no longer satisfies Equation (11). This is logical, since the means of the $\beta_t$'s are, with model M₂, functions of time rather than set to zero, allowing better fitting to the data. As a consequence, the temporal ratios $\lambda_W$ and $\lambda_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₁: $\rho_{\beta^a\beta^W}$ and $\rho_{\beta^a\beta^S}$ 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 ($\beta_t^a, \beta_t^W, \beta_t^S$) time series from model M₂ graphically distinguishable from those of model M₁ displayed in Figure 4.

The values obtained for $\mu_a, \mu_W$ and $\mu_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).

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** PDF of residuals for linear model M₁ in annual (a) and winter (b) point mass-balance series.</figcaption></figure>
The years of change found for winter and summer balances are significantly different. The posterior distribution of $\tau_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.

$\sigma_{\beta w_1}$ is only slightly lower than $\sigma_{\beta w_2}$, so $\tau_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 $\sigma_{\beta w_2}$ is much higher than $\sigma_{\beta w_1}$. $\tau_s$ is therefore also a change point in variance separating two periods with interannual variabilities differing from the overall variance $\sigma_{\beta s}$ 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 $\beta_t$'s 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 ($\tau_w, \tau_s$).

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 $\tau_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 $\beta_t^w$ series of any temporal structure.

For winter balance, the structured trend is strongly discontinuous around $\tau_w = 1976$, indicating that the change is sudden. Between 1949 and 1976, the mean glacier-wide winter budget $\frac{1}{\tau_w - t_0 + 1} \sum_{t_0}^{\tau_w} (\beta_t^w + \langle B^w \rangle)$ is 1.48 m w.e. $a^{-1}$, but is not uniformly distributed since there is a decrease of -0.37 m w.e. $a^{-1}$ over this period ($b_{w1} = -0.013$ m w.e. $a^{-2}$). In 1977, winter mass input increases strongly and its mean $\frac{1}{t_0 + T - \tau_w - 1} \sum_{\tau_w + 1}^{t_0 + T - 1} (\beta_t^w + \langle B^w \rangle)$ between 1977 and 2007 is 1.88 m w.e. $a^{-1}$. Note though that, as observed before 1976, winter mass balance slowly decreases in this second period by 0.20 m w.e. $a^{-1}$ until 2007 ($b_{w2} = -0.007$ m w.e. $a^{-2}$). Mean glacier-wide winter balance was therefore lower than the long-term mean $\langle B^w \rangle = 1.69$ m w.e. $a^{-1}$ 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.

<table><caption>**Table 3.** Posterior characteristics, model M₂ applied on Glacier de Sarennes point mass-balance series for the 1949–2007 period.</caption><thead><tr><th></th><th>Mean</th><th>SD</th><th>2.5%</th><th>97.5%</th></tr></thead><tbody><tr><td>⟨B<sup>3</sup>⟩ (m w.e. a<sup>-1</sup>)</td><td>-0.97</td><td>0.020</td><td>-1.01</td><td>-0.93</td></tr><tr><td>⟨B<sup>5</sup>⟩ (m w.e. a<sup>-1</sup>)</td><td>1.69</td><td>0.014</td><td>1.67</td><td>1.72</td></tr><tr><td>⟨B<sup>8</sup>⟩ (m w.e. a<sup>-1</sup>)</td><td>-2.66</td><td>0.019</td><td>-2.70</td><td>-2.62</td></tr><tr><td>σ<sub>a</sub> (m w.e. a<sup>-1</sup>)</td><td>0.30</td><td>0.015</td><td>0.273</td><td>0.331</td></tr><tr><td>σ<sub>w</sub> (m w.e. a<sup>-1</sup>)</td><td>0.214</td><td>0.011</td><td>0.195</td><td>0.236</td></tr><tr><td>σ<sub>s</sub> (m w.e. a<sup>-1</sup>)</td><td>0.431</td><td>0.021</td><td>0.391</td><td>0.476</td></tr><tr><td>ρ<sub>b</sub><sup>w</sup></td><td>0.39</td><td>0.06</td><td>0.27</td><td>0.50</td></tr><tr><td>ρ<sub>b</sub><sup>w</sup></td><td>0.50</td><td>0.05</td><td>0.42</td><td>0.60</td></tr><tr><td>ρ<sub>b</sub><sup>w</sup></td><td>-0.23</td><td>0.03</td><td>-0.30</td><td>-0.16</td></tr><tr><td>σ<sub>βw1</sub> (m w.e. a<sup>-1</sup>)</td><td>0.464</td><td>0.076</td><td>0.339</td><td>0.636</td></tr><tr><td>σ<sub>βw2</sub> (m w.e. a<sup>-1</sup>)</td><td>0.403</td><td>0.069</td><td>0.296</td><td>0.565</td></tr><tr><td>σ<sub>β1</sub> (m w.e. a<sup>-1</sup>)</td><td>0.608</td><td>0.085</td><td>0.466</td><td>0.796</td></tr><tr><td>σ<sub>β2</sub> (m w.e. a<sup>-1</sup>)</td><td>0.906</td><td>0.145</td><td>0.673</td><td>1.241</td></tr><tr><td>σ<sub>β</sub> (m w.e. a<sup>-1</sup>)</td><td>1.044</td><td>0.020</td><td>1.006</td><td>1.083</td></tr><tr><td>σ<sub>β</sub> (m w.e. a<sup>-1</sup>)</td><td>0.485</td><td>0.015</td><td>0.456</td><td>0.513</td></tr><tr><td>σ<sub>β</sub> (m w.e. a<sup>-1</sup>)</td><td>0.893</td><td>0.019</td><td>0.856</td><td>0.929</td></tr><tr><td>ρ<sub>β</sub><sup>w</sup></td><td>0.52</td><td>0.03</td><td>0.47</td><td>0.57</td></tr><tr><td>ρ<sub>β</sub><sup>w</sup></td><td>0.89</td><td>0.01</td><td>0.87</td><td>0.90</td></tr><tr><td>ρ<sub>β</sub><sup>w</sup></td><td>0.07</td><td>0.04</td><td>0.00</td><td>0.14</td></tr><tr><td>λ<sub>a</sub></td><td>0.92</td><td>0.01</td><td>0.91</td><td>0.94</td></tr><tr><td>λ<sub>w</sub></td><td>0.84</td><td>0.02</td><td>0.80</td><td>0.87</td></tr><tr><td>λ<sub>s</sub></td><td>0.81</td><td>0.02</td><td>0.78</td><td>0.84</td></tr><tr><td>μ<sub>a</sub></td><td>0.17</td><td>0.06</td><td>0.06</td><td>0.30</td></tr><tr><td>μ<sub>w</sub></td><td>0.19</td><td>0.08</td><td>0.05</td><td>0.35</td></tr><tr><td>μ<sub>s</sub></td><td>0.27</td><td>0.08</td><td>0.12</td><td>0.42</td></tr><tr><td>τ<sub>w</sub> (year)</td><td>1976</td><td>2.0</td><td>1974</td><td>1983</td></tr><tr><td>τ<sub>s</sub> (year)</td><td>1982</td><td>2.2</td><td>1977</td><td>1985</td></tr><tr><td>a<sub>w1</sub> (m w.e. a<sup>-1</sup>)</td><td>-0.032</td><td>0.169</td><td>-0.372</td><td>0.293</td></tr><tr><td>a<sub>w2</sub> (m w.e. a<sup>-1</sup>)</td><td>0.494</td><td>0.284</td><td>-0.100</td><td>1.015</td></tr><tr><td>a<sub>s1</sub> (m w.e. a<sup>-1</sup>)</td><td>0.014</td><td>0.185</td><td>-0.355</td><td>0.377</td></tr><tr><td>a<sub>s2</sub> (m w.e. a<sup>-1</sup>)</td><td>0.107</td><td>0.311</td><td>-0.501</td><td>0.719</td></tr><tr><td>b<sub>w1</sub> (m w.e. a<sup>-2</sup>)</td><td>-0.013</td><td>0.012</td><td>-0.034</td><td>0.013</td></tr><tr><td>b<sub>w2</sub> (m w.e. a<sup>-2</sup>)</td><td>-0.007</td><td>0.006</td><td>-0.019</td><td>0.006</td></tr><tr><td>b<sub>s1</sub> (m w.e. a<sup>-2</sup>)</td><td>0.024</td><td>0.011</td><td>0.004</td><td>0.046</td></tr><tr><td>b<sub>s2</sub> (m w.e. a<sup>-2</sup>)</td><td>-0.015</td><td>0.008</td><td>-0.030</td><td>0.000</td></tr></tbody></table>

Regarding summer balance, the transition lasts ~4–5 years around $\tau_s = 1982$. During these 4–5 years, it is shifted down by about 1.11 m w.e. $a^{-1}$. Before 1982, mean glacier-wide summer balance is $\frac{1}{\tau_s - t_0 + 1} \sum_{t_0}^{\tau_s} (\beta_t^s + \langle B^s \rangle) = -2.26$ m w.e. $a^{-1}$, that is 0.4 m w.e. $a^{-1}$ more than the long-term mean $\langle B^s \rangle$. Furthermore, during this first period, summer glacier-wide balance increases by 0.62 m w.e. $a^{-1}$ ($b_{s1} = +0.024$ m w.e. $a^{-2}$). After 1982, summer balances are significantly lower, $\frac{1}{t_0 + T - \tau_s - 1} \sum_{\tau_s + 1}^{t_0 + T - 1} (\beta_t^s + \langle B^s \rangle) = -3.22$ m w.e. $a^{-1}$, i.e. nearly 0.6 m w.e. $a^{-1}$ 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^{-2}$), with summer glacier-wide balance more in deficit by 0.30 m w.e. $a^{-1}$ between $\tau_s$ and 2007.

For glacier-wide annual balance, three regimes are separated by the two change points, $\tau_w$ and $\tau_s$:

1. Before $\tau_w = 1976$, Glacier de Sarennes loses mass with a mean annual wide balance of -0.80 m w.e. $a^{-1}$, but the trend goes from strongly negative balances to less negative balances from year to year ($b_a = +0.011$ m w.e. $a^{-2}$), mainly due to less negative summer balances.

2. Between 1976 and 1982, glacier-wide winter balances increase strongly while the summer regime has not yet
<figure><img src="image_8.png" /><figcaption>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).</figcaption></figure>

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. After the summer balance shift down in $\tau_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. DISCUSSION

## 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 Sareennes series. We have tested a simpler model, M₃, with a single change point, $\tau$, in the two series. The main mode in the distribution of $\tau$ is then in 1982 (Fig. 7c), i.e. the year obtained for $\tau_s$ with model M₂ (Fig. 7a). This year is identified since it is the most significant of the ($\tau_w$, $\tau_s$) pair. The change-point selection is, however, less clear than with model M₂ because $\tau = 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 Sareennes data.

We have also analysed the effect of setting a change point $\tau_a$ in the annual balance instead of the change point $\tau_s$ in the summer balance. This results in a well-marked main mode in 1984 in $\tau_a$'s posterior PDF, thus slightly later than $\tau_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 $\tau_w \neq \tau_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 $\gamma_i^a \delta_t^a$, for the annual balance and $\gamma_i^w \delta_t^w$ for the winter balance in model M₂, leading to the full bivariate ANOVA model M₅:

$$ \begin{pmatrix} b_{it}^a \\ b_{it}^w \end{pmatrix} \sim N_2 \left( \left( \begin{pmatrix} \alpha_i^a + \beta_i^a + \gamma_i^a \delta_t^a \\ \alpha_i^w + \beta_i^w + \gamma_i^w \delta_t^w \end{pmatrix}, \Sigma \right) \right). \quad (33) $$

Among all $\gamma_i$'s, only $\gamma_i^w$ is significantly nonzero at the 95% credibility level for Sareennes. 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 ($\delta_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 $\beta_t$'s (Fig. 9). This confirms that most of the temporal structure is actually captured by the linear temporal terms. Note that the $\delta_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
<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

time is more adapted to the data for annual balance than for winter balance.

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; Spiegelhalter 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 $\sigma_{a\text{cross}}^2 = \frac{1}{NT} \sum_{i=1}^N \sum_{t=t_0}^{t_0+T-1} (\gamma_i^a \delta_i^a)^2$ 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 $\beta_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, Sarenes 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 (Elsberg and others, 2001). However, this effect has no impact here because the $\beta_t^a$ 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 (Vincent, 2002; Huss and others, 2008). The rise in winter precipitation since 1977 that has been found has a more regional signification.
<figure><img src="image_10.png" /><figcaption>Fig. 9. Box plots of the additional nonlinear temporal terms in the annual (a) and winter (b) balances, model M5. The posterior mean is shown in black, the interquartile range in grey and the 95% credibility interval by dashed lines. Temporal cross terms δ_t^a and δ_t^w are displayed in physical units, whereas spatial terms, γ, are arbitrarily normalized to unity.</figcaption></figure>

Indeed, no specific trend is reported in glacier winter balances in the eastern Swiss Alps, for instance on Silvretta and Clariden glaciers (Müller and Kappenberger, 1991; Huss and Bauder, 2009). Nevertheless, closer to Sarennes, Aletschgletscher, Valais, Switzerland, displays a synchro- nous change point in winter balance around 1974–76, though less significant in magnitude than for Sarennes (Huss and Bauder, 2009). Furthermore, this change point is present in various snow and weather series in the French Alps (Durand and others, 2009a,b) and avalanche runout eleva- tion series (Eckert 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 (Thibert and Vincent, 2009), a com- bination 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 M2 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 M2 (before τ<sub>w</sub> = 1976, between τ<sub>w</sub> and τ<sub>s</sub> = 1982, and after τ<sub>t</sub>), 1981 is approximately the middle of the 6–7 year long 1977–84 period of high mass-balance regime.

As indicated by Thibert and Vincent (2009), because of the linearity of the variance decomposition, the glacier-wide balance, B_t^a, can be evaluated from the annual deviation, β_t^a, and the mean annual balance, Φ, obtained from photo- grammetry according to Figure 2b:

$$B_t^a = \Phi + \beta_t^a. \qquad (34)$$

The upper-case letters (B_t^a, Φ) 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 α_t^a terms. The cumulative balance, ΣB_t^a, for glaciological year t is then given by:

$$\Sigma B_t^a = \sum_{q=t_0}^{t} B_q^a = (t - t_0 + 1)\Phi + \sum_{q=t_0}^{t} \beta_q^a. \qquad (35)$$

Under model M2, β_t^a is centred on a linear time trend with slope b_a, which is homogeneous to a rate. The sum of the β_t^a
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Cumulative glacier-wide balance, $\Sigma B_t^a$, obtained by combining photogrammetry performed in 1952, 1981 and 2003, and annual $\beta_t^a$ terms centred on a linear time trend given by model M2. As a result of the two linear trends in annual balances, trends in cumulative balance, $\Sigma B_t^a$, are parabolic over the 1949–81 and 1981–2007 periods.</figcaption></figure>

terms then has a second-order time dependence such as:

$$ \sum_{q}^{t} (a_{a.} + b_{a.}q) = a_{a.}t + \frac{1}{2}b_{a.}t^2 + C, \quad (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:

$$ \frac{\partial^2 \Sigma B_t^a}{\partial t^2} = \frac{\partial B_t^a}{\partial t} = \frac{\partial a^*(t)}{\partial t} = b_{a.} = b_{w.} + b_{s.} \quad (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 M2 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.

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.

## 6. CONCLUSION AND OUTLOOKS

In this paper, a hierarchical modelling framework has been developed to study the spatio-temporal variability of the three components of a glacier mass-balance series, so as to take into account the different seasonal processes. Bayesian estimation was used because of its convenience for hierarchical modelling, but without recourse to informative priors, so as to extract an unperturbed temporal signal from the data. With regard to previous more empirical modelling procedures, this constitutes an important step forward in the study of mass-balance series since it allows an appropriate extraction of refined temporal patterns.

Another important feature of our approach is that it uses a complete bivariate dataset of measured surface balances as input, so as to capture the mean effect at the glacier scale. The extracted temporal signal is therefore free of any effect related to the geometry (surface) and the dynamics of the glacier that affect surface integrated balances, except the retroaction due to elevation changes discussed in section 5.2. It therefore offers an advantageous alternative to the arbitrary selection of a single site, which may discard a part of the climatic information present in the overall data and also lead to misleading conclusions because of local effects.

In detail, a fully multivariate extension was first proposed for Lliboutry's well-known linear ANOVA model. Even if the assumption of Gaussian residuals remains questionable for winter balance, application to the exceptional Sarennes data series has shown that a linear variance decomposition satisfactorily represents the different components of a glacier mass balance. Furthermore, the explicit quantification of all variability sources and correlations was found useful to properly represent the uncertainty around each estimate, and to quantify missing values based on information transferred between sites/years/variables.

Next, different non-exchangeable structures capable of capturing change points and various types of systematic trends were considered. Hierarchy was shown to take into account the uncertainty related to the $\beta_t$'s for their estimation. Different change points in 1976 (highly discontinuous) and 1982 (much more transitional) separating different
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Summary of glacier-wide annual balance $B_t^a = \langle B^a \rangle + \beta_t^a$ decomposition using hierarchy. Empirical estimates are plotted against annual estimates and mean trend from model M$_2$, and respective 95% credibility intervals.</figcaption></figure>

mean trends for respectively winter balance and summer balance as well as different variances for summer balance were highlighted. They explain the complex patterns existing in the corresponding annual mass-balance series.

This hierarchical model is able to identify predominant structured patterns, in our case linear trends with break points, among other interannual fluctuations modelled as Gaussian white noises. Thus, ~20% of the interannual variability of Sarenes winter balance series is related to two linear time trends separated by the change point in 1976, while the rest of the temporal variability is captured by a constant white noise before and after the change point. A relatively analogous feature is identified for summer balance, with 27% of variance related to two linear trends separated by the 1982 change point, but with a white noise which is significantly higher over more recent years. The fraction of the signal captured by the linear trends is not that high, but neither is it negligible. It is presumably of the correct order of magnitude for a marked linear drift due to climate change (including the retroaction due to elevation change) with regard to the truly random interannual fluctuations and to other systematic effects such as periodic fluctuations, long memory effects, multiple change points, etc., which cannot be captured with our model. Furthermore, it is sufficient to well represent the predominant quasi-parabolic pattern in the cumulative balance, for instance the accelerated mass loss over the last 25 years.

As a final illustration of the value of our hierarchical model for decomposing the data in two steps, Figure 11 compares, for the glacier-wide balance, the annual empirical estimates $\frac{1}{N}\sum_{i=1}^{N} b_{it}^a$, to the two components provided by model M$_2$: the annual estimates $\beta_t^a + \langle B^a \rangle$, and the mean trend, $a^*(t) + \langle B^a \rangle$, as well as their 95% credibility intervals. The interannual variability of empirical estimates is clearly very well captured by the annual estimates provided by the model, with a low smoothing effect due to the capacity of the hierarchical model to transfer information from one year to another. The extraction of the mean trend can be considered as a second smoothing that distinguishes the

predominant pattern from the interannual variability, i.e. the distinct behaviours before $\tau_w$ and after $\tau_s$ separated by ~6 years of high mass-balance regime.

Although primarily methodological, our model construction approach was guided by the structure of the studied data. For instance, two change points in the seasonal components of the balance were found to be much better than only one in fitting the data variability, and parameterization in terms of winter and summer balance for the time trends gave much more realistic results than a direct analysis of the actually observed couple, winter and annual balances. Finally, even if the interaction terms were able to capture a small part of the data variability, they were shown to be relatively unnecessary, especially for the study of the $\beta_t$ series which is the primary purpose of this paper. Nevertheless, things could be different for other case studies for which the different proposed model structures M$_1$–M$_5$ could be compared and tested.

Further work could test different modelling extensions. For the temporal trend, more advanced time-series models including nonlinear trends and/or multiple change points instead of only one could be employed, so as to be able to distinguish other types of systematic changes from the 'true' random fluctuations. However, it must be kept in mind that a compromise between model complexity and robustness must be found. Indeed, application to the Sarenes data has shown that, even with our relatively simple model, we are close to the significance limit for several parameters because of the difficulty of extracting information from hidden variables. It would therefore appear to be difficult to expand the model much without losing confidence in the obtained patterns, thus making their climatic relevance questionable.

Other modelling improvements could involve explicit modelling of the spatial structure between the $\alpha_i$'s using geostatistical methods and/or specific elevation effects. An additional hierarchical level could also be considered in the variance decomposition, so as to simultaneously process data from several glaciers instead of only one as previously investigated by Letréguilly and Reynaud (1990) with
Lliboutry's model. These latter points could lead to a more accurate distinction between true spatial gradients and orographic effects, as well as between inter- and intra-glacier spatial variabilities.

However, starting from the proposed basis, the priority is now to bridge the gap with the analysis of direct climatic data as mentioned in section 5.2 in order to reinforce the interpretation of our results in terms of links with climate change. Indeed, since our approach is purely statistical and time-explicit, it has, when used alone, nearly no predictive power. Evaluating the chances of persistence of the highlighted trends beyond 2007 so as to extrapolate the evolution of the different seasonal components of the balance beyond the last studied data therefore involves introducing additional information and/or physical constraints into the analysis. This could be done by using seasonal climatic covariates in the proposed statistical model and/or using its results to constrain and calibrate a physical model of glacier evolution.

## ACKNOWLEDGEMENTS

We are grateful to the many people who have contributed to the constitution of the exceptional data series for Glacier de Sarennes, especially L. de Crécy and F. Valla. Thank to their efforts, data collection has been, despite limited funds, nearly exhaustive for 60 years. The last years of the record were funded by GLACIOCLIM (http://www-lgge.ujf-grenoble.fr/ServiceObs/) through Observatoire des Sciences de l'Univers de Grenoble (OSUG) and Institut des Sciences de l'Univers (INSU). The authors also thank R. Hock and two anonymous reviewers whose comments led to an improved paper.

## REFERENCES

Banerjee, S., B.P. Carlin and A.E. Gelfand. 2004. *Hierarchical modeling and analysis for spatial data*. Boca Raton, FL, Chapman & Hall.

Berger, J.O. 1985. *Statistical decision theory and Bayesian analysis. Second edition*. New York, Springer-Verlag.

Bloomfield, P. 1992. *Trends in global temperature*. *Climatic Change*, **21**(1), 1–16.

Booth, N.B. and A.F.M. Smith. 1982. A Bayesian approach to retrospective identification of change-points. *J. Econ.*, **19**(1), 7–22.

Brooks, S. 1998. Markov chain Monte Carlo method and its application. *J. R. Stat. Soc., Ser. D*, **47**(1), 69–100.

Brooks, S.P. and A. Gelman. 1998. General methods for monitoring convergence of iterative simulations. *J. Comput. Graph. Stat.*, **7**(4), 434–455.

Carlin, B.P., A.E. Gelfand and A.F.M. Smith. 1992. *Hierarchical Bayesian analysis of changepoint problems*. *Appl. Stat.*, **41**(2), 389–405.

Clark, J.S. 2005. *Why environmental scientists are becoming Bayesians*. *Ecol. Lett.*, **8**(1), 2–14.

Clark, J.S. and A. Gelfand, eds. 2006. *Applications of computational statistics in the environmental sciences: hierarchical Bayes and MCMC methods*. Oxford, etc., Oxford University Press.

Diaz, J. 1982. Bayesian detection of a change of scale parameter in sequences of independent gamma random variables. *J. Econ.*, **19**(1), 23–29.

Durand, Y., M. Laternser, G. Giraud, P. Etchevers, B. Lesaffre and L. Mérindol. 2009a. Reanalysis of 44 years of climate in the French Alps (1958–2002): methodology, model validation,

climatology, and trends for air temperature and precipitation. *J. Appl. Meteorol. Climatol.*, **48**(3), 429–449.

Durand, Y., G. Giraud, M. Laternser, P. Etchevers, L. Mérindol and B. Lesaffre. 2009b. Reanalysis of 47 years of climate in the French Alps (1958–2005): climatology and trends for snow cover. *J. Appl. Meteorol. Climatol.*, **48**(12), 2487–2512.

Dyurgerov, M.B. and J.D. Dwyer. 2001. The steepening of glacier mass balance gradients with northern hemisphere warming. *Z. Gletscherkd. Glazialgeol.*, **36**(1–2), 107–118.

Eckert, N., E. Parent, L. Bélanger and S. Garcia. 2007. *Hierarchical Bayesian modelling for spatial analysis of the number of avalanche occurrences at the scale of the township*. *Cold Reg. Sci. Technol.*, **50**(1–3), 97–112.

Eckert, N., E. Parent, M. Naaim and D. Richard. 2008. Bayesian stochastic modelling for avalanche predetermination: from a general system framework to return period computations. *Stoch. Environ. Res. Risk Assess.*, **22**(2), 185–206.

Eckert, N., H. Baya and M. Deschatres. 2010a. Assessing the response of snow avalanche runout altitudes to climate fluctuations using hierarchical modeling: application to 61 winters of data in France. *J. Climate*, **23**(12), 3157–3180.

Eckert, N., M. Naaim and E. Parent. 2010b. Long-term avalanche hazard assessment with a Bayesian depth-averaged propagation model. *J. Glaciol.*, **56**(198), 563–586.

Eckert, N., E. Parent, R. Kies and H. Baya. 2010c. A spatio-temporal modelling framework for assessing the fluctuations of avalanche occurrence resulting from climate change: application to 60 years of data in the northern French Alps. *Climatic Change*, **101**(3–4), 515–553.

Elsberg, D.H., W.D. Harrison, K.A. Echelmeyer and R.M. Krimmel. 2001. Quantifying the effects of climate and surface change on glacier mass balance. *J. Glaciol.*, **47**(159), 649–658.

Funk, M., R. Morelli and W. Stahel. 1997. Mass balance of Griesgletscher 1961–1994: different methods of determination. *Z. Gletscherkd. Glazialgeol.*, **33**(1), 41–55.

Gaucherel, C., F. Campillo, L. Misson, J. Guiot and J.-J. Boreux. 2008. Parameterization of a process-based tree-growth model: comparison of optimization, MCMC and Particle Filtering algorithms. *Environ. Model. Softw.*, **23**(10–11), 1280–1288.

Geman, S. and D. Geman. 1984. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. *IEEE Trans. Patt. Anal. Mach. Intell.*, **6**(6), 721–741.

Gilks, W.R., S. Richardson and D.J. Spiegelhalter, eds. 1996. *Markov chain Monte Carlo in practice*. London, Chapman & Hall.

Haeberli, W. 1995. *Glacier fluctuations and climate change detection*. *Geogr. Fís. Din. Quat.*, **18**(2), 191–199.

Hoinkes, H. 1970. Methoden und Möglichkeiten von Massen-haushaltsstudien auf Gletschern: Ergebnisse der Messreihe Hintereisferner (Ötztaler Alpen) 1953–1968. *Z. Gletscherkd. Glazialgeol.*, **6**(1–2), 37–90.

Huss, M. and A. Bauder. 2009. 20th-century climate change inferred from four long-term point observations of seasonal mass balance. *Ann. Glaciol.*, **50**(50), 207–214.

Huss, M., A. Bauder, M. Funk and R. Hock. 2008. Determination of the seasonal mass balance of four Alpine glaciers since 1865. *J. Geophys. Res.*, **113**(F1), F01015. (10.1029/2007JF000803.)

Kaser, G., J.G. Cogley, M.B. Dyurgerov, M.F. Meier and A. Ohmura. 2006. Mass balance of glaciers and ice caps: consensus estimates for 1961–2004. *Geophys. Res. Lett.*, **33**(19), L19501. (10.1029/2006GL027511.)

Kuhn, M. 1984. Mass budget imbalances as criterion for a climatic classification of glaciers. *Geogr. Ann.*, **66**(A3), 229–238.

Lauritzen, S.L. 1996. *Graphical models*. Oxford, etc., Oxford University Press.

Letréguilly, A. and L. Reynaud. 1990. Space and time distribution of glacier mass-balance in the Northern Hemisphere. *Arct. Alp. Res.*, **22**(1), 43–50.

Lliboutry, L. 1974. Multivariate statistical analysis of glacier annual balances. *J. Glaciol.*, **13**(69), 371–392.
Mayo, L.R., M.F. Meier and W.V. Tangborn. 1972. A system to combine stratigraphic and annual mass-balance systems: a contribution to the International Hydrological Decade. *J. Glaciol.*, **11**(61), 3–14.

Mearns, L.O., C. Rosenzweig and R. Goldberg. 1997. Mean and variance change in climate scenarios: methods, agricultural applications, and measures of uncertainty. *Climatic Change*, **35**(4), 367–396.

Meier, M.F. and W.V. Tangborn. 1965. Net budget and flow of South Cascade Glacier, Washington. *J. Glaciol.*, **5**(41), 547–566.

Müller, H. and G. Kappenberger. 1991. Claridenfirn-Messungen 1914–1984. *Zürcher Geogr. Schr.* 40.

Nogaj, M., S. Parey and D. Dacunha-Castelle. 2007. Non-stationary extreme models and a climatic application. *Nonlinear Proc. Geophys.*, **14**(3), 305–316.

Oerlemans, J. and J.P.F. Fortuin. 1992. Sensitivity of glaciers and small ice caps to greenhouse warming. *Science*, **258**(5079), 115–117.

Oerlemans, J. and N.C. Hoogendoorn. 1989. Mass-balance gradients and climatic change. *J. Glaciol.*, **35**(121), 399–405.

Ohmura, A., A. Bauder, H. Müller and G. Kappenberger. 2007. Long-term change of mass balance and the role of radiation. *Ann. Glaciol.*, **46**, 367–374.

Parent, E. and J. Bernier. 2007. *Le raisonnement bayésien: modélisation et inférence*. Paris, etc., Springer-Verlag.

Perreault, L., J. Bernier, B. Bobée and E. Parent. 2000a. Bayesian change-point analysis in hydrometeorological time series. Part 1. The normal model revisited. *J. Hydrol.*, **235**(3–4), 221–241.

Perreault, L., J. Bernier, B. Bobée and E. Parent. 2000b. Bayesian change-point analysis in hydrometeorological time series. Part 2. Comparison of change-point models and forecasting. *J. Hydrol.*, **235**(3–4), 242–263.

Rao, A.R. and W. Tirtotjondro. 1996. Investigation of changes in characteristics of hydrological time series by Bayesian methods. *Stoch. Hydrol. Hydraul.*, **10**(4), 295–317.

Rasmussen, L.A. 2004. Altitude variation of glacier mass balance in Scandinavia. *Geophys. Res. Lett.*, **31**(13), L13401. (10.1029/2004GL020273.)

Rasmussen, L.A. and L.M. Andreassen. 2005. Seasonal mass-balance gradients in Norway. *J. Glaciol.*, **51**(175), 601–606.

Solomon, S. and 7 others, eds. 2007. *Climate change 2007: the physical science basis. Contribution of Working Group I to the*

Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge, etc., Cambridge University Press.

Soruco, A. and 9 others. 2009. Mass balance of Glaciar Zongo, Bolivia, between 1956 and 2006, using glaciological, hydrological and geodetic methods. *Ann. Glaciol.*, **50**(50), 1–8.

Spiegelhalter, D.J., N.G. Best, B.P. Carlin and A. van der Linde. 2002. Bayesian measures of model complexity and fit. *J. R. Stat. Soc.*, **64**(4), 583–639.

Tanner, M.A. 1992. *Tools for statistical inference: observed data and data augmentation methods*. Berlin, etc., Springer.

Thibert, E. and C. Vincent. 2009. Best possible estimation of mass balance combining glaciological and geodetic methods. *Ann. Glaciol.*, **50**(50), 112–118.

Thibert, E., R. Blanc, C. Vincent and N. Eckert. 2008. Glaciological and volumetric mass-balance measurements: error analysis over 51 years for Glacier de Sarennes, French Alps. *J. Glaciol.*, **54**(186), 522–532.

Valla, F. 1989. Forty years of mass-balance observations on Glacier de Sarennes, French Alps. *Ann. Glaciol.*, **13**, 269–272.

Valla, F. and C. Piedallu. 1997. Volumetric variations of the Glacier de Sarennes, French Alps, during the last two centuries. *Ann. Glaciol.*, **24**, 361–366.

Vallon, M. and J.C. Leiva. 1981. Bilans de masse et fluctuations récentes du Glacier de Saint-Sorlin (Alpes Françaises). *Z. Gletscherkd. Glazialgeol.*, **17**(2), 143–167.

Vallon, M., C. Vincent and L. Reynaud. 1998. Altitudinal gradient of mass-balance sensitivity to climatic change from 18 years of observations on glacier d'Argentière, France. *J. Glaciol.*, **44**(146), 93–96.

Vincent, C. 2002. Influence of climate change over the 20th century on four French glacier mass balances. *J. Geophys. Res.*, **107**(D19), 4375. (10.1029/2001JD000832.)

Vincent, C. and M. Vallon. 1997. Meteorological controls on glacier mass balance: empirical relations suggested by measurements on glacier de Sarennes, France. *J. Glaciol.*, **43**(143), 131–137.

Vincent, C., G. Kappenberger, F. Valla, A. Bauder, M. Funk and E. Le Meur. 2004. Ice ablation as evidence of climate change in the Alps over the 20th century. *J. Geophys. Res.*, **109**(D10), D10104. (10.1029/2003JD003857.)

Von Storch, H. and F.W. Zwiers. 2002. *Statistical analysis in climate research*. Cambridge, Cambridge University Press.

Wikle, C.K. 2003. Hierarchical Bayesian models for predicting the spread of ecological processes. *Ecology*, **84**(6), 1382–1394.

MS received 10 June 2010 and accepted in revised form 15 October 2010