1. Introduction
The spatial distribution of specific balances of a glacier measured at individual sites is a key parameter used to optimize the sampling site network and obtain the glacier-total balance. Spatial dependence of mass balance has been studied in various papers mostly in relation to altitude dependence (Reference LliboutryMeier and Tangborn, 1965; Reference Gilks, Richardson and SpiegelhalterHoinkes, 1970; Khun, 1984; Reference Oerlemans and HoogendoornRasmussen, 2004; Reference RasmussenRasmussen and Andreassen, 2005). These authors all concluded that altitude profiles, b(z, t), of balances measured at year t vary from year to year by an amount Δb(t) that is constant over the entire glacier. This can be written:
where α(z) is a specific function for each glacier that can be integrated over the altitude range of the glacier to give the glacier-total balance. However, some studies also report a deviation from Equation (1) (Reference Meier and TangbornOerlemans and Hoogendoorn, 1989; Reference Fountain and VecchiaFunk and others 1997; Dyurgerov and Dwyer, 2001). The linear decomposition of temporal and spatial variability of Equation (1) has been generalized by Reference KuhnLliboutry (1974) from measurements over the ablation area of Glacier de Saint-Sorlin. Reference VallaVallon and Leiva (1981) extended later the relation to the accumulation area. They proposed the following equation:
whee2re b(x, y, t) is the mass balance recorded at location (x, y) for year t and the α term represents the spatial effects at (x, y). The α(x, y) term can be reduced to α(z) when the altitude variable alone can be used to describe the spatial dependence of balances. The β(t) term represents the systematic balance differences from year to year over the entire glacier and is analogous to the amount Δb(t) over which the vertical profile is shifted annually in Equation (1). For a given period of record of mass balances measured on a sampling network, the mean balance measured at (x,y) must be represented by one of the two components of the variance decomposition of Equation (2). Under the constraint ββ(t) = 0, the mean balance over the whole period of record is given by α(x, y) while β(t) represents the systematic annual differences from the mean and is defined as the centred balance.
We present here a complete analysis of the mass-balance series of Glacier de Sarennes, France, to check if the annual balance exhibits a linear variance decomposition, as previously demonstrated for a part of the period of record (Reference Thibert, Faure and VincentThibert and others, 2008). Another objective is to use the variability analysis results to determine the best possible estimation of the glacier-total balance using an independent measurement of the cumulative mass balance obtained from photogrammetry (geodetic mass balance).
2. Data
Mass-balance series has been measured on Glacier de Sarennes (45°07’ N; 6°07’ E), Massif des Grandes Rousses, France, since 1949 (Reference Thibert, Blanc, Vincent and EckertValla, 1989). Sarennes is a small glacier (0.4km2 in 2003) with a small altitude range of approximately 200m (2840–3050ma.s.l. in 2003). The annual balance 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. Balances are measured at sites 1–5 for which the mean elevations over the period of record are, respectively, 2860, 2920, 2945, 2995 and 3020ma.s.l. (Fig. 1). 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–2007 data table (59 years × 5 sites = 295 data) is therefore an incomplete dataset (27 missing values; Table 1).
Fig. 1. Location of stakes where balances are measured at Glacier de Sarennes. Complementary sites where balances have been measured occasionally on the sides of the glacier are also indicated. Contour lines on the glacier are at 25 m intervals.
Table 1. Table of mass balances, bi,t , measured at the five sites, centred balance, βt , deduced from the variance analysis, and glacier-total balance, bt , obtained by combining the geodetic and glaciological methods. Boldface values are estimated missing values from the variance analysis
3. Variance Analysis
A general variance decomposition analysis (CETAMA, 1986) can be used to determine whether the mass balance follows Equation (2). The idea is to separate spatial and temporal effects according to the mathematical model (Reference KuhnLliboutry, 1974):
where bi , t is the mass balance recorded at stake 〉 for year t.
The α〉 term represents the spatial effects at location i, and their mean corresponds to the mean annual mass balance at this location over the whole period, while βt represents the annual deviation from the mean (ββt = 0). The cross-term, βiβt , accounts for non-linear effects that deviate from the time and space variable separation under the constraints ββt = 0 and ββ〉 = 0. If cross-terms are negligible, Equation (3) is reduced to a pure linear model (α〉 and βt terms). The єi ,t term represents residuals corresponding to both measurement errors and discrepancies between the model and data (unexplained variance). The temporal control factors and residuals are assumed to be identically normally distributed with means equal to zero. The spatial (σα ), temporal (σβ ), coupled (σβ ) and residual (σє ) variance components are to be estimated. As the experience table is incomplete, the missing values are treated as additional variables to be estimated. Inference is carried out in a Bayesian framework using a Markov chain Monte Carlo scheme (Reference Funk, Morelli and StahelGilks and others, 1996).
4. Results
The variance analysis results are presented in Table 2. The variance explained by the decomposition according to Equation (3) is high (96%), showing that the adopted variance analysis model is well suited for Glacier de Sarennes. It makes it possible to examine separately spatial, temporal and non-linear effects, as well as residuals. Whether the modelling hypotheses are fulfilled is checked by standard tests such as Snedecor tests (variance comparison) and 2 tests (normality of distributions) applied to each component of the decomposition as described in the following subsections.
Table 2. Variance analysis results over the period 1949–2007. The explained variance is the variability (%) explained by the decomposition of variance according the two possible models (with cross-terms and linear). Taking into account non-linear effects would only explain an additional variance of 1%
4.1. Spatial variability
The spatial variability of the data is σα = 0.25mw.e. a–1. The deviation from the mean caused by sampling at five sites is significant, showing that sampling on more than one site is worthwhile. However, sampling on four sites would be acceptable, as the discrepancy between the α〉 spatial terms is not significant for stakes 1 and 3 (Table 3). As the spatial sampling at Sarennes is only located on the central axe of the glacier (Fig. 1), whether such an adopted sampling covers the overall spatial variability of the mass balance over the glacier has been investigated in different ways:
Table 3. Spatial effects at each sampling site (αi ), standard deviation of residuals (σє ), stake-to-βt correlations (ri 2 ) and mean weighting coefficients (si ) used in the conventional area integration relationship. Best estimator in the estimation of βt is shown in bold. The mean altitude of each stake over the period of record is also indicated
1. First, some complementary mass balances have been measured over 3–6 year long periods on the sides of the glacier (Fig. 1). Such measurements do not show significant discrepancy with stakes 3 and 4 (within the measurement errors). These additional data could have been included in the variance analysis to check if they increase the spatial variability. However because these additional data cover a very short time period, estimation of missing values would have led to variance analysis results of too high uncertainty to give any firm conclusion on that point.
2. Second, as explained in a previous study (Reference Thibert, Faure and VincentThibert and others, 2008), the static glacier hypothesis (i.e. there is negligible ice flow) can be adopted at Sarennes, so that the water equivalent of local altitude changes can be considered as mass balances. Altitude changes obtained from the 1952–2003 photogrammetry do not show spatial variation on the sides of the glacier in comparison to the central axis.
3. Third, according to Reference Thibert, Faure and VincentThibert and others (2008), the spatial variability of the mass balance deduced from photo-grammetry and the static glacier hypothesis is 0.27mw.e. This value is very close to the 0.25mw.e. obtained from the variance analysis of the five axial sites.
These three points suggest that the adopted sampling covers the overall spatial variability of the mass balance at Sarennes.
4.2. Temporal variability
The temporal variability of the annual balance (σβ = 1.05mw.e. a–1) is the main source of variability in the data (90%). The centred balance has a normal distribution (Fig. 2). As Equation (3) fits the data well, the mass balance at each stake differs from year to year by an amount that is nearly uniform over the glacier. This is illustrated in Figure 3 where the mean balance at each stake (α〉 curve) has been plotted versus altitude. For a given year of the period of record (e.g. 1995; Fig. 3), the altitudinal profile is the same as the mean profile of α〉 but shifted from the annual deviation β 1995 = +1.67mw.e. a–1. This linear behaviour is not surprising for Sarennes, which is a small glacier with a limited altitude range, and such a linear law has already been reported over the period 1952–2003 (Reference Thibert, Faure and VincentThibert and others, 2008), as well as for other glaciers (Reference LliboutryMeier and Tangborn, 1965; Reference KuhnLliboutry, 1974; Reference HoinkesKuhn 1984; Reference Oerlemans and HoogendoornRasmussen, 2004).
Fig. 2. Annual deviation of the mass balance (centred balance, βt ) and its temporal variability (σβ term). βt is the annual amount (constant all over the glacier) from which altitudinal profiles are shifted from year to year in Figure 3. The hypothesis of normal distribution of βt terms is acceptable.
Fig. 3. Mean annual balances at each site (α〉 spatial terms) and their spatial variability (σα term). The linear composition of variance is illustrated with the altitudinal profile recorded in 1995, which is the same profile as the mean α〉 profile but shifted by an amount of β 1995 = +1.67mw.e.
4.3. Cross-terms
Non-linearity accounts for only 1% of the variance of the data (σβ = 0.12mw.e. a–1). Figure 4 shows cross-terms estimated for the period 1958–2002 where there is no missing value in the data. Fifteen values are significantly different from zero (at 95% confidence level). There is a small spatial effect in non-linearity, as cross-terms have a significantly higher variance at stake 1 (σβ = 0.24mw.e. a–1). This can be explained by topographic shading and enhanced winter accumulation due to snowdrift in this area. Regarding the complete period 1949–2007, the estimation of missing values in the data (first 8 years of the period of record; Table 1) increases noise in the analysis: confidence intervals of β cross-terms increase from 0.46 to 0.78mw.e. This means that cross-terms are not significantly different from zero. Moreover, the variance of cross-terms (0.12mw.e.a–1) is low compared with random errors of measurement of mass balances (0.17mw.e.a–1 at Sarennes; Reference Thibert, Faure and VincentThibert and others, 2008). Considering their low variance, non-linear effects are neglected hereafter and the decomposition is considered as purely linear (95% of explained variance).
Fig. 4. Non-linear terms (γiςt ) as a function of time and their variability (12 cmw.e. a–1) over the period 1958–2002. The dotted line shows the cross-terms at stake 1, where the variability is higher.
4.4. Residuals
The residuals (at each site and overall) have normal and centred distributions (Fig. 5). The variance of residual is indicated in Table 3. Because differences of residuals distributions are not significant between different sites, residuals are not correlated with space. Similarly, residuals are not correlated temporally (only 5.6% of common variance), as Figure 6 shows, even if a trend with a slight increase is detectable over the period of record. For the model to be fully acceptable, the residuals of uncontrolled factors must be compared with measurement errors (0.17mw.e. a–1 at Sarennes). The standard deviation of the residuals (Table 2) is therefore just slightly greater than the uncertainty of field measurements. As the difference between the residual variance and the measurement variance can be considered as a rough estimate of the model error, this indicates that the model fits the data well.
Fig. 5. Residuals of the linear model at each stake. Stake 4 is the site that displays the lowest residuals (σє = 19 cmw.e.) and which is the best estimator of the centred balance (see Fig. 7). Highest residuals are observed at site 1 (σє = 34 cmw.e.) which also displays the highest non-linear terms (see Fig. 4).
Fig. 6. Overall residuals of the linear model (including the variance of the non-linear terms) as a function of time. There is a trend with slight increase over the period of record.
Considering negligible cross-terms, the temporal and spatial decomposition that can be adopted for the annual balance t at stake 〉 is given by a pure linear model that explains 95% of the variance of the data. This results in the following numerical terms:
where the residual 
accounts for the additional variance (1%) of the neglected cross-terms.
5. Glacier-Total Balance
The spatial integration of balances measured at individual sites can be done using classical regression methods with respect to altitude (Reference CogleyCogley, 1999; Reference Dyurgerov and DwyerFountain and Vecchia, 1999) when balances are correlated with altitudes, or using geostatistical tools such as the kriging model or Thiessen polygons (Reference VallaVallon and Leiva, 1981; Reference Fountain and VecchiaFunk and others, 1997) when variograms show spatial correlation of balances. At Sarennes, for any year t in the N-year period of record, the glacier-total balance is calculated from a linear combination of point measurements weighted by the surface ratio, si , associated with each site i:
where the s〉 terms are calculated from the mean area distribution with altitude over the period of record with the constraint βs〉 = 1 (Table 3). As bi , t follows the decomposition of spatial and temporal variables of Equation (3), Equation (5) gives:
where the sum Ps〉 α〉 corresponds to the mean annual balance integrated over the glacier area which can be defined as ʃ α, and where the sum βs〉 βt equals βt because βt terms are functions of the time variable t alone and βs〉 = 1. As a combination of centred Gaussian distributions, the last term of Equation (6) is itself Gaussian with an expected value of 0. Equation (6) is mainly important in that the annual deviation at each site equals the annual deviation in the overall bt , so that estimating the annual glacier-total balance can be performed in two steps: (a) estimating βt terms each year; and (b) estimating the constant ʃ α which is specific to the glacier.
5.1. Estimation of the βt terms
Because βt is an annual uniform response of all the stakes as shown by Equations (2) and (3), only one stake is needed to record the annual βt deviation. Best estimators for this are stakes that show smallest residuals, єi , t in Equation (3), and the best correlation with βt . Stake-to-βt correlations, ri , are shown in Table 3. Stake-to-βt correlations are effectively high (at least >89.9%). The best estimator is stake 4 (and to a lesser extent stake 2) which explains 97.4% of the annual response of all stakes (Fig. 7).
Fig. 7. Estimation of the centred balance, βt , from a single stake. Stake 4 is the best estimator of the centred balance, as it explains alone 97.4% of the variance of the centred balance.
5.2. Estimation of the ʃ α term
Over the period of record of N years, the cumulative balance is given by:
Because ʃ α is constant with time, the sum βʃ α is N ʃ α and the sum ββt is 0 under the constraints of the linear model (section 3). Therefore the α integration term can be expressed as:
As a mean of a centred population, the last term β єi , t /N of Equation (8) is equal to 0, and therefore an unbiased estimator of ʃ α is the mean annual value of the glacier-total cumulative balance over the period of record. A consequence of the separation of temporal and spatial variables is therefore the possibility of using an independent estimation, BN/N, of the spatial integration of the α〉 terms over the glacier area when calculating the annual glacier-total balance. An alternative to Equation (6) is therefore to substitute ʃ α and write:
where BN is the balance obtained from photogrammetry over the same period (geodetic mass balance). Because it is based on bedrock fixed reference control points, while glaciological mass balances are referred to the previous year late summer surface, and because it is measured over a longer time period, the cumulative balance obtained from photogrammetry is less affected by systematic error. Moreover, photogrammetry integrates the overall spatial variability of the mass balance over the glacier surface and consequently is a better estimator than the integration of the α〉 terms (Reference Thibert, Faure and VincentThibert and others, 2008).
Regarding random errors, the variance of b phot t obtained from Equation (9) is for a fixed year t:
where σBN is the error associated with the balance obtained from photogrammetry, and αs〉 the (unknown) error associated with the weighting coefficients of the area integration.
The variance of bt obtained from Equation (6) is:
so that the variances of the two possible estimations of the glacier-total balance only differ by the first right-hand terms of Equations (10) and (11). Error associated with the estimation of the cumulative mass balance from photogrammetry is typically 1 m w.e. (Reference Rasmussen and AndreassenThibert and others, 2005, 2008)which gives 
with N = 59 years in Equation (10). Using σα
= 0.25m w.e. a–1 in Equation (11), the square root of the sum 
is 13 cm w.e. A 10% error in the weighting coefficients would result in an additional contribution, 
and gives a close value of random error of 14 cmw.e. a–1 in Equation (11). Calculation according to Equation (6) seems therefore as much affected by random errors as Equation (9).
Thus, as an equivalent variance and unbiased estimator of the mean annual value of the glacier-total cumulative balance, preference is given to photogrammetry to estimate the spatial component in the glacier-total balance according to:
where the random error term on the right-hand side, ±0.19 mw.e., is just slightly above the annual mean error of ±0.17 mw.e. obtained with conventional area integration (Reference Thibert, Faure and VincentThibert and others, 2008) and accounts for random errors in photogrammetry (through the σBN term) and in field data (through the unexplained variance σє term).
The glacier-total balance that we recommend from the data is therefore a combination of the cumulative balance (BN = –31.5±1.07mw.e.) measured from photogrammetry for the period 1952–2003 (N = 51; Reference Thibert, Faure and VincentThibert and others, 2008), and the centred balance obtained from the variance analysis (Table 1). As no volumetric balance is available for the initial and final periods of the record (8 years among 59 years), we have used the glacier-total balance obtained from the glaciological mass balance (Equation (6)).
6. Conclusion
The purpose of this study was to extend and confirm the variance analysis developed by Reference Thibert, Faure and VincentThibert and others (2008) to the entire period of record of the mass-balance series of Glacier de Sarennes.
A first consequence of the separation of spatial and temporal variability is that the mass balance differs from year to year by an amount that is uniform over the glacier, leading to the possibility of recording this annual deviation (centred balance) from a single site to which it is highly correlated (97%). A second consequence of the linear separation of time and space variables is the possibility of estimating the spatial component in the glacier-total balance by an independent method, such as photogrammetry which is less affected by systematic errors than the area integration of balances at each site. However, it should not be concluded that the number of sites used to measure the annual balance can be reduced automatically, because random errors related to sampling and to point measurements are reduced as the number of sites increases and the risk of losing some data is higher for a low number of sites. Moreover, for larger glaciers with a higher altitude range, time variability is expected to be lower at higher altitudes (Reference Vallon and LeivaVallon and others, 1998) so that stake-to-centred balance correlations are expected to be lower. High correlations for Sarennes are likely related to the small size of the glacier. Such a method combining glaciological and geodetic methods could be one of the standard approaches while the glaciological community progresses towards homogeneous mass-balance series.
Acknowledgements
We thank L. de Crécy and F. Valla who have maintained the long mass-balance records of Glacier de Sarennes since 1949. We thank N. Eckert for help with the variance analysis. We are grateful to two anonymous reviewers who helped to improve the quality of the paper.
