Reanalysing glacier mass balance measurement series

Abstract. Glacier-wide mass balance has been measured for more than sixty years and is widely used as an indicator of climate change and to assess the glacier contribution to runoff and sea level rise. Until recently, comprehensive uncertainty assessments have rarely been carried out and mass balance data have often been applied using rough error estimation or without consideration of errors. In this study, we propose a framework for reanalysing glacier mass balance series that includes conceptual and statistical toolsets for assessment of random and systematic errors, as well as for validation and calibration (if necessary) of the glaciological with the geodetic balance results. We demonstrate the usefulness and limitations of the proposed scheme, drawing on an analysis that comprises over 50 recording periods for a dozen glaciers, and we make recommendations to investigators and users of glacier mass balance data. Reanalysing glacier mass balance series needs to become a standard procedure for every monitoring programme to improve data quality, including reliable uncertainty estimates.

where a is the magnitude of the misalignment, b is the direction of the misalignment, and c is the mean elevation bias (z adjustment ) divided by the tangent of the mean slope of the sample. The input data (dh, α and ω) should be a sample from non-glacierized stable terrain that contains a distribution of at least half of all possible aspects, uniformly distributed, if possible. Terrain slope and aspect are ideally computed from the master DEM, however, since the process is iterative, the choice of DEM for slope and aspect is arbitrary. The unknown co-registration parameters (a, b, and c) can be solved using a least-squares minimization, common to programmes such as Excel, MATLAB and R. Since this is an analytical solution based on a non-analytical terrain surface, the process should be repeated until the solution converges. Thus, once the initial solution is determined, the slave DEM should be translated by the magnitude (a) and direction (b) of the co-registration vector and adjusted for the mean vertical bias, and the DEMs re-differenced. The translation of the parameters , , and into x, y, and z adjustments is as follows: Equation A1 can then be solved again. Typically no more than two to three iterations are required depending upon the quality of the terrain sample.
If 3D co-registration is successful, then the bias of the co-registration (z adjustment in Eq. A4) is removed and there remains the uncertainty of the vertical co-registration adjustment(s). One approach to estimate this potential un-removed vertical error is to introduce additional elevation datasets, either as a control or as a part of a time series. When three or more datasets are available, co-registration can be performed between each of these, and the summation of the 3D co-registration vectors returns the residual error remaining within the series. Studies have shown that remaining vertical errors can reach magnitudes of at least 1-3 m (Nuth et al., 2012;Berthier et al., 2012) and should be included in the uncertainty assessment.
There are a few exceptions in which the solution of the 3D co-registration problem fails. The first is on flat or low-sloped terrain; i.e., slopes less than three to five degrees that are present in the stable terrain sample. Second, many older maps and DEMs created for glaciological applications do not contain a sufficient sample of the surrounding topography. In these cases, co-registration may be performed using alternate (image matching) methods, described in Berthier et al. (2007), for example.

Supplement B.
A method to determine the spatial auto-correlation based on Rolstad et al. (2009).
The uncertainty in the spatially averaged elevation difference is estimated by following these steps: 1. Create an elevation difference grid of the bedrock region surrounding the glacier. 2. Detrend the grid, if necessary, using a polynomial model to remove bias (as described in Supplement A). 3. Estimate the spatial auto-correlation, such as by statistically assessing the grid to determine the semi-variogram parameters of nugget c 0 , partial sill c 1 and range a 1 by fitting a spherical semi-variogram model to the empirically derived semi-variogram. Standard geostatistical software packages are available to do this. For reference, the standard deviation of the elevation error derived over bedrock σ ∆z is related to the semi-variogram parameters by: 4. If the correlation range a 1 is greater than the representative radius L of the averaging area where Δh is the pixel size.
If the correlation range is less than the representative radius of the averaging area (a 1 <L), as may be the case when determining the geodetic mass balance over large areas, then σ S is determined using 2 = 0 ∆ℎ 2 2 + 1 5 1 1 2 2 a 1 < L .
As discussed in Rolstad et al. (2009) there may be more than one scale of spatial correlation related to the derivation of the DEMs. It must be emphasized that it is generally the largest correlation scale that has the greatest impact on the spatially averaged uncertainty.