Best possible estimation of mass balance combining glaciological and geodetic methods

Abstract

Variance analysis of the long time series of mass balances recorded on Glacier de Sarennes (45°07’ N, 6°07’ E) France, since 1949 shows that the variability can be separated linearly in two spatial and temporal terms. Annual balances deviate from their mean values over the period of record by an annual amount that is uniform over the glacier. Annual balances at each stake are therefore highly correlated, and sampling at a single site would be acceptable to record the annual deviation. A result of the linear character of the variance is the possibility of obtaining a systematic error-free estimate of the annual glacier-total budget by combining the mean annual balance obtained from photogrammetry and the annual deviation obtained from the variance analysis, rather than using the traditional area integration of balances at each stake.

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 (Meier and Tangborn, 1965; Hoinkes, 1970; Khun, 1984; Rasmussen, 2004; Rasmussen 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:

(1)

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) (Oerlemans and Hoogendoorn, 1989; Funk and others 1997; Dyurgerov and Dwyer, 2001). The linear decomposition of temporal and spatial variability of Equation (1) has been generalized by Lliboutry (1974) from measurements over the ablation area of Glacier de Saint-Sorlin. Vallon and Leiva (1981) extended later the relation to the accumulation area. They proposed the following equation:

(2)

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 (Thibert 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 (Valla, 1989). Sarennes is a small glacier (0.4km² 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, b_i,t , measured at the five sites, centred balance, β_t , deduced from the variance analysis, and glacier-total balance, b_t , 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 (Lliboutry, 1974):

where b_i , _t is the mass balance recorded at stake 〉 for year t.

(3)

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 (Gilks 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 ² 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 (r_i ² ) and mean weighting coefficients (s_i ) 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 (Thibert 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 Thibert 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 β ₁₉₉₅ = +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 (Thibert and others, 2008), as well as for other glaciers (Meier and Tangborn, 1965; Lliboutry, 1974; Kuhn 1984; Rasmussen, 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 β ₁₉₉₅ = +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; Thibert 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:

(4)

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 (Cogley, 1999; Fountain and Vecchia, 1999) when balances are correlated with altitudes, or using geostatistical tools such as the kriging model or Thiessen polygons (Vallon and Leiva, 1981; Funk 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, s_i , associated with each site i:

(5)

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 b_{i , t} follows the decomposition of spatial and temporal variables of Equation (3), Equation (5) gives:

(6)

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 b_t , 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, r_i , 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:

(7)

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:

(8)

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, B_N/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:

(9)

where B_N 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 (Thibert and others, 2008).

Regarding random errors, the variance of b ^phot _t obtained from Equation (9) is for a fixed year t:

(10)

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 b_t obtained from Equation (6) is:

(11)

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. (Thibert 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:

(12)

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 (Thibert 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 (B_N = –31.5±1.07mw.e.) measured from photogrammetry for the period 1952–2003 (N = 51; Thibert 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 Thibert 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 (Vallon 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.