An 11-year record of mass balance of Brewster Glacier, New Zealand, determined using a geostatistical approach

Abstract

Recognising the scarcity of glacier mass-balance data in the Southern Hemisphere, a mass-balance measurement programme was started at Brewster Glacier in the Southern Alps of New Zealand in 2004. Evolution of the measurement regime over the 11 years of data recorded means there are differences in the spatial density of data obtained. To ensure the temporal integrity of the dataset a new geostatistical approach is developed to calculate mass balance. Spatial co-variance between elevation and snow depth allows a digital elevation model to be used in a co-kriging approach to develop a snow depth index (SDI). By capturing the observed spatial variability in snow depth, the SDI is a more reliable predictor than elevation and is used to adjust each year of measurements consistently despite variability in sampling spatial density. The SDI also resolves the spatial structure of summer balance better than elevation. Co-kriging is used again to spatially interpolate a derived mean summer balance index using SDI as a co-variate, which yields a spatial predictor for summer balance. The average glacier-wide surface winter, summer and annual balances over the period 2005–15 are 2484, −2586 and −102 mm w.e., respectively, with changes in summer balance explaining most of the variability in annual balance.

1. INTRODUCTION

Monitoring surface mass balance of mountain glaciers is important as fluctuations in glacier mass over time have hydrological implications for water resource users and are used to infer changes in climate. An understanding of the contribution of glaciers and ice caps is also crucial in any assessment of global sea-level rise. A globally coordinated effort to obtain mass-balance measurements has enabled glacier-wide results from >250 glaciers to be documented (Zemp and others, 2015), with >4000 annual observations available spanning the past seven decades. However, the dataset is strongly biased towards the Northern Hemisphere and Europe. There are only 30 glaciers that have (almost) continuous measurements since 1976 (Zemp and others, 2009) and of these only one (Echaurren Norte, Andes) is located in the Southern Hemisphere, which is small (area < 0.40 km²) compared with the average size of the present 30 reference glaciers (7.25 km² in 2003).

To address the issue of Southern Hemisphere glaciers being heavily under-represented in the global glacier monitoring network, a mass-balance programme on Brewster Glacier in the Southern Alps of New Zealand was initiated in 2004. The glaciers in the Southern Alps are unique by virtue of their location in the south-west Pacific, situated mid-way between Antarctica and the Equator. The Southern Alps are surrounded by vast areas of ocean and are strongly influenced by both subtropical and polar air masses, with the interaction of these contrasting air masses in the prevailing westerly airflow resulting in synoptic scale atmospheric circulation having a strong influence on glacier behaviour (Fitzharris and others, 1992, 1997; Clare and others, 2002; Gillett and Cullen, 2011). The mass-balance programme on Brewster Glacier has provided a platform from which to assess glacier sensitivity to atmospheric forcing in the Southern Alps (Anderson and others, 2010; Gillett and Cullen, 2011; Conway and Cullen, 2013, 2016; Conway and others, 2015; Cullen and Conway, 2015).

Despite the importance of the glaciological record from Brewster Glacier, there have been numerous challenges associated with the task of producing a high-quality time series of mass balance. The challenges we have encountered on Brewster Glacier are no different to those reported by others contributing to the global monitoring network, with the main sources of random and systematic error linked to uncertainties associated with field measurements at point locations, spatial averaging of these observations over an entire glacier and the treatment of changes in area and elevation over time (e.g. Zemp and others, 2013). To ensure that the glaciological observations and a high quality mass-balance record from Brewster Glacier are accessible to the glaciological community, the main purpose of this research is to describe a new geostatistical method developed to produce an 11-year record of mass balance. We anticipate our new approach will provide a methodological framework that allows historical ‘debates’ such as the relevancy of the number of point measurements taken each year in a mass-balance programme (e.g. Fountain and Vecchia, 1999) to be scrutinised in a more systematic way.

We hypothesise that the spatial structure of winter and summer balance on Brewster Glacier for any individual mass-balance year contains sufficient similarity to other years to allow spatially distributed indices to be developed for each seasonal balance. To test this hypothesis we use a new geo-statistical approach to develop unique indices to calculate winter and summer mass balance. We also use the approach to create a predictor for annual balance, which allows us to revisit years when summer balance was not obtained. The glacier-wide estimates of winter, summer and annual mass balance are used to identify changes in equilibrium line altitude (ELA) and accumulation area ratio (AAR), as well as to resolve mass-balance gradients over the observational period. Finally, we assess the relationship between the mass-balance record and a longer record of ELA determined from the end-of-summer snowline (EOSS) programme (Chinn and others, 2012), which has been the benchmark for monitoring glaciers in the Southern Alps.

2. SITE DESCRIPTION

Brewster Glacier is a small, maritime glacier situated in the Southern Alps immediately west of the main divide (44.08°S, 169.43°E), with a surface area of 2.03 km² (Fig. 1). The length of the glacier is ~2.5 km, with its elevation ranging between 1700 and 2400 m a.s.l. At higher altitudes (>2000 m a.s.l.), the glacier has a south-westerly aspect and a mean slope of 31°, while the lower, larger portion of the glacier (77% of its total area) has a southerly aspect and is more gently sloping (mean slope of 10°). The upper and smaller part of the glacier situated below the summit of Mount Brewster contains a number of large ice cliffs. Brewster Glacier is one of 50 ‘index glaciers’ in the Southern Alps that has been surveyed as part of the EOSS programme (Chinn and others, 2005, 2012; Willsman and others, 2015). The first oblique aerial photography obtained for Brewster Glacier as part of the EOSS programme was captured in 1978, and it now provides a 38-year record of glacier behaviour. Prior to 1999, end-of-summer snowlines were generally lower than the long term average (positive mass balance), while higher snowlines (negative mass balance) have been observed more frequently since 2008, a trend consistent with other index glaciers in the Southern Alps.

Fig. 1. The ablation stake network, probing and snow pit locations on Brewster Glacier. The glacier margins are shown for 1986 and 2011. The 2011 margin was digitized from a Quickbird image from 8 February 2011 and late summer aerial photographs on 13 and 30 March 2011. It should be noted that the probing locations shown above 2000 m a.s.l. are not regularly sampled, while those below typically are (see text for further details). Coordinates shown are defined using the New Zealand Transverse Mercator 2000 projection (NZTM2000).

3. GLACIER OUTLINE AND SURFACE TOPOGRAPHY

Glacier outlines and the surface topography of Brewster Glacier have been produced and refined on several occasions. The DEM most commonly used to characterise the surface topography of Brewster Glacier is based on a GPS survey of the glacier tongue and parts of the accumulation area carried out in January 1997 (Willis and others, 2009). Surface topography for the upper accumulation area was obtained by digitising 38 additional points from a Land Information New Zealand (LINZ) 1:50 000 topographic map (NZMS260/G38), which was based on images obtained in a photogrammetric aerial survey in 1986. Willis and others (2009) used the same topographic map to define the glacier boundary, which has subsequently been updated using GPS measurements of the terminus in conjunction with the mass-balance programme. Surface topography data obtained from GPS measurements are available from surveys conducted in 1997, 2005 and 2012 but the products do not contain full coverage of Brewster Glacier.

For the purpose of this study, a 15 m resolution DEM interpolated from 20-m contours derived from the photogrammetric aerial survey in 1986 (NZSoSDEM, Columbus and others, 2011) was used for all calculations involving elevation due to its complete representation of the surface. Comparison between the 15-m resolution DEM and the 2012 GPS survey exhibited elevation differences ranging from −30 m at the terminus to +30 m in the accumulation area. The glacier boundary considered representative of the observation period was defined using a near snow-free orthorectified (1.6 m RMSE) pan-sharpened Quickbird image from the KiwImage project acquired on 8 February 2011. To aid in the identification of the glacier boundary, oblique aerial photographs captured on 13 and 30 March 2011 that contained very little seasonal snow were also used. By using a fixed geometry, the glacier-wide mass-balance record we present is a reference-surface mass balance (Elsberg and others, 2001). The area of the glacier decreased from 2.54 km² in 1986 to 2.03 km² in 2011, equivalent to a 21% decrease over the 25-year period (Fig. 1).

4. GLACIOLOGICAL DATASET

The current mass-balance programme on Brewster Glacier was initiated in February 2004, with the first year of measurements described by George (2005). As the programme evolved changes were made to the sampling methodology, resulting in variability in the glaciological dataset. To clarify how this has been managed, the main characteristics of the dataset are described in the following two sections. Further details about the programme over the period 2004–08 can be found in Stumm (2011), while Figure 1 shows the ablation stake and snow probing locations for the present mass-balance network. The winter mass balance referred to in this research represents the sum of accumulation and ablation over the winter season (Cogley and others, 2011), with the expectation that a small proportion of winter accumulation is lost to ablation. Likewise, the summer mass balance refers to the sum of accumulation and ablation over the summer season, with the ablation of snow and ice being the dominant component. The time system in which mass balance is calculated is from field surveys on floating calendar dates (a floating-date system), thus the duration of the mass-balance year varies slightly from year to year (Table 1).

Table 1. Observational data for each mass balance (MB) year, including the observation type, start and end date of measurements, the method used to obtain the data and the number of measurements. Each MB year contains three rows of information, with the density and depth data used primarily for winter balance, while the stake data are either used for summer or annual balance

^a Annual balance calculated using the method described by Stumm (2011; see pg. 284–285).

4.1. Winter mass-balance data

At the beginning of the mass-balance programme, winter balance was determined at the end of the primary accumulation period (early to mid-November) from three snow pits along the central flowline of the glacier between 1700 and 1900 m a.s.l. (Fig. 1) and a limited number of snow depth measurements (Stumm, 2011) (Table 1). Snow pits were dug to the previous summer surface and density measurements were obtained. On some occasions snow density measurements were also obtained from a snow pit dug near the terminus. In November 2006, a Kovacs ice corer was also used to obtain an estimate of winter balance in the upper part of the glacier at an elevation of 2300 m a.s.l. (Fig. 1). From November 2007 and onwards an emphasis was placed on determining the spatial variability in snow depth, with extensive geo-located probing of snow cover carried out over the glacier (Fig. 1; probe locations), with at least one snow pit dug to obtain density measurements along the central flowline. Variability in the winter-balance measurements stems from the fact that the number of snow density measurements from snow pits varies from year to year, as does the number of snow depth measurements (Table 1). A key data deficiency is the under-sampling of snow density and depth on the steeper, south-west facing slopes in the upper accumulation area.

4.2. Summer mass-balance data

Summer mass balance has been obtained by measuring changes in the height of a network of PVC stakes drilled into the glacier using a Heucke steam ice drill. On most occasions, ablation stakes were drilled in conjunction with the end of winter field measurements (early to mid-November), and re-surveyed at the end of summer (mid to end of March). The network was initially extensive, with over 30 ablation stakes drilled in December 2004 at 100 m intervals along the central flowline, as well as four transects extending across the central flowline (George, 2005; Stumm, 2011). In November 2008, the stake network was reduced to permit a smaller number of people to complete the field measurements, with further minor changes to their locations made in November 2014 (Fig. 1).

Like the winter mass-balance measurements, there is also sampling variability in the ablation stake data over the observational period (Table 1). The most significant differences in data collection occurred during the 2006–08 mass-balance period, when the network of stakes was used primarily for the calculation of annual balance rather than summer balance (Stumm, 2011). At other times, ablation stakes were not steam drilled at all locations due to time constraints and/or not re-surveyed at the end of summer because they had either melted out of the ice or were completely lost.

5. METHODS

5.1. Calculation of winter mass balance

5.1.1. Snow density

The density of snow at the end of winter was obtained from snow pits from at least one, and up to four, locations on the glacier (Table 1). To assess whether elevation is a control on snow density, a comparison of mean snow densities from all individual snow pits during the programme was made. A one-way analysis of variance (ANOVA) for all years except 2006 and 2010 yields a P-value = 0.26, enabling us to reject the null hypothesis that mean snow densities are significantly different from each other. Only snow densities measured to calculate mass balance in 2007 (P-value < 0.01) and 2011 (P-value < 0.05) display significant differences to other years. The mean snow density from all years is 492 kg m⁻³, with a standard deviation of 36 kg m⁻³. In view of the difficulty of obtaining significant spatial or temporal variations of snow density based on the observational dataset, it was concluded that the best approach to estimate snow density (ρ _s) was to rely solely on the mean snow density $\mu _t^{\rho _s} $ in any given year t and to use an uncertainty of 36 kg m⁻³ for error propagation in the mass-balance calculation.

5.1.2. Creating a snow-depth index

To develop a robust method to estimate winter balance, we searched for a linear predictor capable of explaining the variance in snow depth obtained from the snow pit and snow probing measurements. In alpine terrain, snow depth is usually found to correlate reasonably well with elevation, which is also the case on Brewster Glacier (Fig. 2a). For each year, the snow depth data exhibit a linear correlation with elevation below 2000 m a.s.l. that is significant at the 95% confidence level for all years except 2004 (P-value = 0.26) and 2005 (P-value = 0.10), when the number of snow depth samples obtained were small. On average, elevation is found to explain 73% of the variance in snow depth below 2000 m a.s.l. The small number of data points above 2000 m a.s.l. suggests there is a decline in the snow depth gradient (Fig. 2a), which compromises the linear regression approach to model snow depth across the entire elevation range of Brewster Glacier.

Fig. 2. (a) The relationship between snow depth and elevation. The linear regressions are computed over an elevation range between 1700 and 2000 m a.s.l. (b) Snow depth data for each year are offset using the regressed value at an elevation of 1850 m a.s.l.

Despite this limitation, Figure 2a shows there is a coherency in the spatial structure of snow depth gradient from year to year, which is used as justification to search more thoroughly for a single spatial predictor against which data from individual years might be regressed. To build a model that removes the variability associated with differences in total accumulation from year to year, the linear regression models of each individual year are used to offset the snow depth measurements based on the mid-elevation value of each regression (i.e. 1850 m a.s.l.). The mid-elevation is the mid-point of the elevation range (1700–2000 m a.s.l.) over which the linear regression models between snow depth and elevation are made. Using this approach we subtract the regression value at 1850 m a.s.l. from the snow depth data for each year, which yields a synthetic dataset that shifts all snow depth information to a common origin (Fig. 2b). The consistency of the scatter plot is explained by the limited variance associated with the gradient (mean μ = 11.4 mm m⁻¹, standard deviation σ = 2.3 mm m⁻¹) and suggests that the snow cover at the end of the winter exhibits a coherent spatial structure and is primarily affected by year to year variability in the magnitude of winter balance.

To characterise and model the spatial structure of snow depth, we assume that the offset snow depth data stem from the realisation of a random function (Wackernagel, 2003), which we define here as a snow depth index (SDI). This provides us with a platform from which to apply a multi-variate geostatistical approach to resolve the spatial variability in snow depth. To guarantee first- and second-order stationarity assumptions are met for variogram computation and subsequent kriging (Wackernagel, 2003, p. 51), the SDI and elevation covariates are de-trended using a two-dimensional (2-D) hyperplane. Figures 3a, b demonstrate that the semi-variograms for SDI and elevation exhibit a marked spatial structure, with the SDI range approaching 500 m, while the range in elevation is almost twice as large (~1000 m). The SDI semi-variogram exhibits a nugget of 0.12 m², which corresponds to the collocated variance associated with data from different dates. If all measurements had been collected on the same date, this nugget would capture the measurement error associated with the small-scale effects of snow depth variability (~0.35 m).

Fig. 3. Experimental semi-variograms for (a) SDI and (b) elevation, along with their variogram models. The nugget of the semi-variogram is the intercept at 0-lag. It corresponds to the collocated variance present in the dataset resulting from short-scale variations not captured by the sampling network, measurement errors, and, in the present case, variance inherent to data collected at different dates. The range of the semi-variogram is the distance at which the sill or plateau is reached and is indicative of the loss of spatial auto-correlation between measurements separated by such distance. (c) Cross semi-variogram associated with both variables with fitted linear model of coregionalization. (d) Ergodic cross-correlogram computed for both +h and –h directions to support a complete picture of spatial covariation (Rossi and others, 1992).

The semi-variogram models are computed to enable a linear model of coregionalization (LMC) to be utilised (Goovaerts, 1998), with the purpose of building a co-kriging system capable of predicting SDI across the entire surface of Brewster Glacier. The cross-semivariance and cross-correlogram shown in Figures 3c, d demonstrate that, in addition to the correlation at a glacier-wide scale (Fig. 2), SDI and elevation display a significant spatial cross-correlation structure. By taking advantage of the 15 m resolution DEM that offers a high number of elevation samples, co-kriging can be utilised to obtain a robust prediction of SDI variability, which is expected to result in smaller error variances than univariate kriging or other interpolation methods (Goovaerts, 1998). We used the Geostatistical Analyst Toolbox in ArcGIS v10.2 to apply co-kriging to the SDI samples (Fig. 4a), using the DEM as a cofactor (Fig. 4b). A smooth surface is achieved by using between 5 and 200 samples within a search neighbourhood of 750 m to accommodate the range of both covariates. The outcome of this approach yields a map of prediction and standard error of the regionalised SDI, referred to as SDI _CK and $\sigma _{CK}^{SDI} $ , respectively (Figs 4c, d).

Fig. 4. (a) Set of SDI samples cokriged with elevation from the DEM in (b). (c) Resulting SDI _CK map and (d) associated kriging standard error $\sigma _{CK}^{SDI} $ .

5.1.3. Using the SDI to determine winter mass balance

There is a robust relationship between SDI _CK and snow depth for each winter mass-balance year on record (Fig. 5). The years with only three–four snow depth measurements (2004–06), as determined from snow pits, exhibit significant determination coefficients (P-value < 0.01). All other years have very high and significant correlations (P-values < 10⁻¹¹), which confirm that SDI _CK is a strong predictor of snow depth. On average, SDI _CK explains 89% of the variance in snow depth, which is a considerable improvement compared with using elevation as a predictor (Fig. 2a). Furthermore, the consistent linear relationship between SDI _CK and snow depth does not break down at higher elevations, which is observed to occur when using elevation as the independent variable. Although   $\sigma _{CK}^{SDI} $ remains high due to the distance to support points (Fig. 4d), it is reduced by benefiting from the secondary information of elevation available in this area.

Fig. 5. Measured snow depth versus SDI. Years are offset for clarity.

Our new method to spatially interpolate snow depth is not dependent on data from individual years but instead uses the SDI _CK predictor. The winter balance years with sparse data for which interpolation might have been compromised and subject to artefacts can now be spatially interpolated in a consistent manner. The spatial interpolation for any individual year exploits data from the entire record, which is important in the interpolation of snow depth at higher elevations, while also useful in spreading measurement errors. The adjustment of mass-balance data in this manner shares commonalities with the approaches used by Lliboutry (1974) and Eckert and others (2011), which are examined further in Section 6.7.

Finally, the co-kriging standard error map provides a means to display the statistical uncertainties associated with snow depth that can be propagated to estimate winter balance (b _w). In practice, in any given year (t) and regardless of the sparsity of snow depth data, an estimate of the snow depth random function $\hat D_{\rm t}$ is obtained via linear regression between observations and SDI _CK as:

(1) $$\hat D_{\rm t} = \max \left( {\left\{ {s_{w_t} \; SDI_{CK} + o_{w_t}, 0} \right\}} \right)$$

where s _{w t} and o _{w t} are the slope and intercept of the regression (w, winter balance). This yields a map of winter balance (b _w) for any given year (t) as:

(2) $$b_{w_t} = \mu _t^{\rho _s} \hat D_t$$

The standard error in winter balance for any given year (t) $\sigma _{b_{w_t}}$ is obtained via error propagation:

(3) $$\sigma_{b_{w_{t}}} = \sqrt {{\hat D}_t \times {0.036}^2 + \mu _t^{\rho _s} \times {\left( {s_{w_t} \; \sigma _{CK}^{SDI}} \right)}^2} $$

From the map and standard error of winter balance, the glacier-wide surface winter balance (B _w) and associated standard error $\sigma _{B_w}$ are computed as arithmetic and quadratic means for any given year (t), respectively:

(4) $$\left\{ {\matrix{ {B_{w_t} = \overline {b_{w_t}}} \cr {\sigma _{B_{w_t}} = \sqrt {\overline {\sigma _{b_{w_t}}^2}}} \cr}} \right.\; $$

5.2. Calculation of summer mass balance

5.2.1. Creating a mean summer balance index

Ablation of snow and ice is the main component of the summer balance (b _s ), with measurements from ablation stakes available to calculate b _s for all years except 2006, 2007 and 2008. Summer balance is calculated from observed changes in ablation stake height in relation to the surface for any given year (t) SH _t ( x _i ) at a location x _i to b _s in w.e. Using the mapped snow depth in the previous year $\hat D_{t - 1}$ it becomes possible to estimate b _s as follows:

(5) $$b_{s_t}\left({{\bi x}_i} \right) = \left\{ \matrix{- \mu _{t - 1}^{\rho _s} {\hat D}_{t - 1}\left({{\bi x}_i} \right) - \rho _i\left( SH_t ({{\bi x}_i}) - {\hat D}_{t - 1}({{\bi x}_i}) \right), \cr if\; \; SH_t ( {{\bi x}_i} ) \gt {\hat D}_{t - 1} ({{\bi x}_i}) \cr \quad {- \mu _{t - 1}^{\rho _s} {\hat D}_{t - 1} ({{\bi x}_i}), \; otherwise} \cr} \right.$$

where ρ _i is the density of ice (900 kg m⁻³). Following the same approach as for snow depth, elevation was initially considered as a predictor for summer balance. Correlations between summer balance and elevation proved relatively weak and sometimes not significant, with an average R ² = 54% (Fig. 6a). Alternatively, substantially better correlations are obtained between summer balance and SDI _CK , with an average R ² = 73% (Fig. 6b). Similar to the probing data, a synthetic dataset combining all summer balance estimates is obtained by deducting, for each year, the value of a linear regression model between SDI _CK and summer balance at SDI _CK  = 400 mm. This corresponds approximately to the 1850 m a.s.l. elevation used to offset the snow depth data. However, unlike the snow depth data obtained from snow pit and probing measurements, ablation stakes are always drilled at the same geographic location. As such, all summer balance measurements cannot be used together in a geostatistical approach because of collocation (i.e. lag = 0). Instead, the offset summer balance estimates are averaged to derive a single mean summer balance index (MSBI) at each stake location. The MSBI takes advantage of the statistical averaging, which yields a reduced dispersion and a strong, significant correlation with SDI _CK (R ² = 85%, P-value < 0.01) (Fig. 6b).

Fig. 6. (a) Correlation between summer ablation and elevation at the stakes. (b) Correlation between summer ablation and SDI _CK . Ablation data for each year are offset using the regressed value at elevation SDI _CK  = 400 mm. Black circles are the average ablation at each stake.

Figures 7a, b show the semi-variogram analysis for detrended MSBI and SDI _CK . Both covariates exhibit a marked spatial structure, with a similar range of ~500 m. Spatial co-variation is significant as shown on the cross semi-variogram and cross-correlogram (Figs 7c, d). Given the strong correlation also shown between both variables (Fig. 6b), a co-kriging approach is applied to spatially interpolate MSBI using SDI _CK as a co-variate. This yields a spatial predictor for summer balance (MSBI _CK ), which supports the interpolation of summer balance in regions that do not contain ablation stakes, thus reducing the standard error $\sigma _{CK}^{MSBI} $ (Fig. 8). All samples within a radius of 2000 m are used to interpolate a smooth MSBI _CK surface. Although this radius exceeds the range of both covariates it is required to enable MSBI to be extrapolated over the whole surface but it comes at the expense of generating relatively large errors far from support points.

Fig. 7. Experimental semi-variograms for (a) MSBI and (b) SDI _CK , along with their variogram models. (c) Cross semi-variogram associated with both variables with fitted linear model of coregionalization. (d) Ergodic cross-correlogram computed for both +h and –h directions to support a complete picture of spatial covariation (Rossi and others, 1992).

Fig. 8. (a) Set of MSBI samples cokriged with SDI _CK in (b). (c) Resulting MSBI _CK map and (d) associated kriging standard error $\sigma _{CK}^{MSBI} $ .

5.2.2. Using MSBI to determine summer mass balance

The regression of summer balance values at ablation stake locations versus MSBI _CK for each year is shown in Figure 9. It demonstrates that the linear relationship between MSBI _CK and summer balance is robust for individual years, with all having strong and significant correlations (P-values < 0.01). This verifies that MSBI _CK is a robust predictor of summer balance, explaining on average 78% of the variance, which is vastly improved compared with only using elevation.

Fig. 9. Summer balance versus MSBI. Years are offset for clarity.

Using the same approach as for snow depth (and winter balance), summer balance is obtained via regression of the ablation stake observations over an individual summer with MSBI _CK . This allows an adjustment of the observations and a distribution of observational errors, while minimizing inter- and extrapolation artefacts. For any given year (t) a map of summer balance $(b_{s_t})$ is obtained via linear regression between observations and MSBI _CK as:

(6) $$b_{s_{t}} = \min \left( {\left\{ {s_{s_t} \; MSBI_{CK} + o_{s_t}, 0} \right\}} \right)$$

where s _{s t} and o _{s t} are the slope and intercept of the regression (s, summer balance). The standard error in summer balance for any given year (t) $\sigma _{b_{s_t}}$ is obtained by:

(7) $$\sigma _{b_{s_{t}}} = s_{s_t} \; \sigma _{CK}^{MSBI} $$

From the map and standard error of summer balance, the glacier-wide surface summer balance (B _s ) and associated standard error $\sigma _{B_s}$ are computed by arithmetic and quadratic mean for any given year (t), respectively:

(8) $$\left\{ {\matrix{ {B_{s_t} = \overline {b_{s_{t}}}} \cr {\sigma _{B_{s_t}} = \sqrt {\overline {\sigma _{b_{s_{t}}}^2}}} \cr}} \right.$$

5.3. Calculation of annual mass balance

5.3.1. Years with summer and winter mass-balance measurements

For all years except 2006, 2007 and 2008 a map of annual balance (b _a) is obtained by summing the map of summer balance (b _s) in year t with that of winter balance (b _w) in year t−1:

(9) $$b_{a_t} = b_{w_{t - 1}} + b_{s_t} $$

The independence of errors in b _w and b _s is tested by examining the residuals of cross validation of the co-kriging models for SDI _CK and MSBI _CK . The sets of cross-validation residuals for both spatial predictors on collocated support points (n = 34) fail the Anderson-Darling test (P-value = 0.98 and 0.49, respectively), indicating that both sets of residuals are from a population with a normal distribution. An insignificant Pearson correlation coefficient (r = −0.17, P-value = 0.34) and no visible pattern in the scatter plot could be taken to suggest independence between both sets of residuals. However, independence between both sets of residuals is confirmed at the 99% confidence level by applying the Hilbert-Schmidt Independence Criterion (HSIC, Gretton and others, 2005) (test statistic = 0.2702 < test threshold = 0.5424). The standard error in the mass-balance map for any given year (t) is therefore obtained by Gaussian error propagation by:

(10) $$\sigma _{b_{a_t}} = \sqrt {\sigma _{b_{s_t}}^2 + \sigma _{b_{w_{t - 1}}}^2} $$

Similarly, the glacier-wide surface mass balance (b _a) for any year (t) is obtained as B _{a t}  = B _{w t−1}  + B _{s t} , with a standard error calculated as:

(11) $$\; \sigma _{B_{a_t}} = \sqrt {\sigma _{B_{s_t}}^2 + \sigma _{B_{w_{t - 1}}}^2} $$

5.3.2. Years with only winter and annual mass-balance measurements

In 2006, 2007 and 2008 the network of ablation stakes was only used to support the calculation of annual mass balance (end of summer to the next), which is described in detail by Stumm (2011). In order to incorporate these years into the time series of mass balance, a predictor for annual balance in the form of an annual balance index (ABI) is computed using:

(12) $$ABI = \rho _s\; SDI_{CK} + MSBI_{CK}$$

with an associated standard error:

(13) $$\sigma _{ABI} = \sqrt {{\left( {\rho _s\; \sigma _{CK}^{SDI}} \right)}^2 + \sigma {_{CK}^{MSBI}} ^2} $$

where ρ _s (snow density) is assumed to be equal to 500 kg m⁻³. The linear relationships between measurements of annual balance and the values calculated using ABI for each individual year between 2006 and 2008 are shown in Figure 10. Strong and significant correlations are obtained for 2006 and 2007 (P-values < 0.01). It is encouraging that mass-balance observations during these 2 years can be determined with a high degree of confidence using ABI, which is derived from measurements from all other years, suggesting there is a coherent spatial structure in winter and summer balance. The linear regression for 2008 is not as robust, with a determination coefficient of R ² = 45% (P-value = 0.45). We interpret this as evidence that the proposed approach may also be capable of detecting inconsistencies in the observational record, which are described by Stumm (2011), but we cannot rule out that secondary processes that are not fully accounted for might be having an effect on mass balance (e.g. Purdie and others, 2015). Importantly, the regression approach based on a predictor of annual balance does benefit from redundancy that helps to mitigate errors associated with the observations, while producing a spatially consistent map of annual balance.

Fig. 10. Annual balance versus ABI. Years are offset for clarity.

For the 3 years discussed, a map of annual balance (b _a) for any given year (t) is obtained via linear regression between observations and ABI:

(14) $$b_{a_{t}} = s_{a_t} \; ABI + o_{a_t} $$

where s _{a t} and o _{a t} are the slope and intercept of the regression (a, annual balance). The standard error in annual balance for any given year (t) $\sigma _{b_{a_t}}$ is obtained by:

(15) $$\sigma _{b_{a_{t}}} = s_{a_t} \; \sigma _{ABI}$$

Finally, this allows maps of summer balance to be calculated by subtracting winter balance from annual balance, while the associated errors are obtained by Gaussian error propagation. The average and associated standard errors of the glacier-wide values are computed as arithmetic and quadratic means, respectively.

5.4. ELA and AAR

Positive values on any given annual balance map define the accumulation area, with the boundary from positive to negative values interpreted as the equilibrium line. The accumulation area ratio for any given year (AAR _t ) is computed as the relative share of the accumulation area versus the area of the glacier. The equilibrium line altitude for any given year $\mu _{ELA_t}$ is defined as the average elevation along the equilibrium line. Since elevation does not uniquely control annual balance in our approach, the variance of elevation along the equilibrium line is also captured by $\sigma _{ELA_t}$ . Both AAR _t and $\mu _{ELA_t}$ are compared with estimates from the EOSS programme (Chinn and others, 2012).

6. RESULTS AND DISCUSSION

6.1. Performance of the co-kriging

Cross validation of the co-kriging of SDI and MSBI results in mean errors of 0.11 mm and 4.5 mm w.e., respectively, which are indicative of negligibly biased models. The RMSE is 441 mm for SDI and 211 mm w.e. for MSBI. The largest uncertainty in SDI is on the northern slopes of the accumulation area, with the maximum RMSE equal to 726 mm. This corresponds ~10% error relative to the average snow depth of 7500 mm at this location. At the same location, the RMSE for MSBI reaches a maximum of 251 mm w.e., or 38% of the average summer balance of −650 mm w.e. Although the determination of the co-kriging standard error is rigorous in a geostatistical sense, the quantification of the uncertainty in the upper accumulation area needs to be further scrutinised. The standard error maps are useful as they highlight the need to obtain more data in the upper accumulation area. In the current framework, such samples would have the immediate effect of refining the predictor of snow depth and summer balance, as well as reducing the co-kriging variance, which would result in lowering the uncertainties associated in calculating glacier-wide winter, summer and annual balances.

Not only is the SDI _CK map useful as a tool to calculate winter balance, it provides detailed information about the spatial variability of snow depth. Although the smoothness of the map cannot resolve the fine spatial structure indicative of the expected heterogeneity in snow distribution (Kerr and others, 2013), locations with accumulation, consistently greater or less than average are identified. In particular, the relatively thin snow cover along the central flowline, despite an increase in elevation, raises the hypothesis that snow is removed either by wind and/or is subject to enhanced ablation. Distributed energy and mass-balance modelling (Mölg and others, 2009) or application of a wind field model (e.g. Winstral and Marks, 2002; Dadic and others, 2010) would help to identify the physical processes, beyond those imposed by the hypsometry of the glacier, responsible for the observed spatial variability.

6.2. Mass balance time series for Brewster Glacier

The time series of mass balance for Brewster Glacier during the period 2005–15 are shown in Table 2. The glacier-wide surface winter (B _w), summer (B _s) and annual balances (B _a) averaged over the observation period are 2484, −2586 and −102 mm w.e., respectively, indicating that the glacier has been losing mass over the period. The variability of B _s is more pronounced than that of B _w, with their standard deviations equal to 813 and 372 mm w.e., respectively. There is no significant correlation between B _s and B _w (R ² = 0.08, P-value = 0.39) and they are statistically independent as determined using the HISC (Gretton and others, 2005) (test statistic = 0.2196 < test threshold = 0.4918). There is a weak (R ² = 0.38) but significant (P-value < 0.05) correlation between B _a and B _w. A much greater coefficient of determination is obtained between B _a and B _s (R ² = 0.87, P-value < 0.01), suggesting that changes in summer balance explain most of the variability in annual balance. This finding is consistent with the findings of Anderson and others (2010), as well as being observed, since the mid-20th century, at a large number of other mountain glaciers (Zemp and others, 2015).

Table 2. Winter, summer and annual mass-balance values for Brewster Glacier. The units for mass-balance values are mm w.e., while mean snow density $(\mu _{t - 1}^{\rho _s} )$ for winter balance is kg m⁻³

The values in bold show the mean snow density and mass-balance values.

^a Quadratic average.

6.3. Comparison of annual ELA and AAR values with EOSS data

The EOSS altitude, which is used as a surrogate for ELA, has been independently estimated from oblique aerial photography for the entire observational period (Chinn and others, 2005, 2012). Comparing values of ELA and AAR calculated from our spatial modelling of mass balance using glaciological data to those estimated from EOSS observations shows a significant correlation between values derived by the two independent methods (P-values < 0.01 for ELA and AAR). The spatial modelling approach succeeds in capturing some of the variability in ELA and AAR as determined from the EOSS programme. However, the large ELA discrepancies in 2008 and 2012 and general departure from the ideal 1:1 ratio line, suggest shortcomings in either one or both methods.

Over the 11-year mass-balance period (2005–15), the snowline from the EOSS survey was digitally georeferenced on four occasions, and not always from rectified photographs. For all other years, the snowline altitude is established by applying an interpolation method, which involves arranging all photographs into a sequence based on increasing area of snow cover (descending order of ELAs). ‘The photograph being interpolated is then carefully compared and inserted into its appropriate place in the sequence’ (Willsman and others, 2015, p. 15). Following this approach, 2008, 2011 and 2012 were assessed as having nearly identical ELAs. Review of these photographs confirms that Brewster Glacier is largely snow free, with only patchy snow cover visible on the north-eastern slopes. This snow versus ice pattern is also largely captured by the spatial modelling approach, with relatively high ELA in 2008 and 2011 but not 2012. Although the glaciological record confirms that the winter balances during the mass-balance years of 2008 and 2011 are similar, as are the summer balances, there are substantial differences between 2011 and 2012 (Table 2). To avoid any uncertainty that may have been introduced from our extrapolation method we examined only collocated stakes and probing locations to assess the differences in mass loss and gain in 2011 and 2012, with mean surface height change at collocated ablation stakes equal to 6845 and 4796 mm, respectively. A comparison of snow depth at the end of the 2010 and 2011 winters, determined from collocated snow probing locations, shows that the snowpack was 28% thinner in 2011, with winter balance 20% lower after taking into account variations in snow density. A comparison of annual balance at the collocated ablation stake and snow probing locations shows that mass loss in 2012 was only 35% of that observed in 2011. This confirms that the 2 years were substantially different in terms of annual balance, despite the glacier being assessed as having similar ELAs using the EOSS survey interpolation method.

Removing 2012 from the analysis improves the correlation of ELA derived by the two methods from R ² = 0.63 to R ² = 0.82, with the latter linear regression not significantly different to the 1:1 ratio line at the 95% confidence level (Fig. 11a). Any uncertainty in the determination of ELA following the two methods also has implications for the calculation of AAR. However, when ELA is >2000 m a.s.l., variations in the two estimates only correspond to small variations in AAR because of the hypsometry of the glacier. For this reason, the linear regression between the two differently derived AAR values is more robust (R ² = 0.75, P-value < 0.01) than the ELA relationship over the entire period (Fig. 11b), and further improves if 2012 is discarded (R ² = 0.79). This suggests that AAR may be a better primary variable to assess mass balance than ELA.

Fig. 11. Comparison between (a) ELA and (b) AAR from the spatial modelling of mass balance and the values determined from the EOSS (Willsman and others, 2015).

6.4. Resolving equilibrium state ELA and AAR

To determine values of ELA and AAR on Brewster Glacier representative of the observation period, a mass-balance map of the glacier in an equilibrium state $(b_{a_{eq}})$ is defined as:

(16) $$b_{a_{eq}} = \overline {b_{a_t}} - \overline {B_{a_t}} $$

where $\overline {b_{a_t}} $ is the spatially averaged map computed from the sequence of annual balance maps and $\overline {B_{a_t}} $ the mean of the glacier-wide annual balances. By definition, the glacier in equilibrium $(b_{a_{eq}})$ has a glacier-wide mass balance equal to zero (balanced-budget), from which representative equilibrium values of ELA and AAR can be determined using the $b_{a_{eq}} = 0$ line (Fig. 12a). The uncertainties in ELA_eq and AAR_eq are obtained using the uncertainty associated with $b_{a_{eq}}$ determined by propagating the mapped standard error of $b_{a_t}$ with that of $\overline {B_{a_t}} $ (Table 2), which results in 1898 ±39 m a.s.l. < ELA_eq = 1923 ± 25 m a.s.l. < 1938 ± 28 m a.s.l. and 0.35 < AAR_eq = 0.40 < 0.47 (Fig. 12b).

Fig. 12. (a) Mass-balance map of the glacier in an equilibrium state $(b_{a_{eq}})$ and (b) associated standard error.

Given the importance of using EOSS as a method to monitor glacier behaviour in the Southern Alps (Chinn and others, 2005, 2012), it is of interest to compare ELA_eq and AAR_eq values to those determined from the EOSS programme. As described by Chinn and others (2012), EOSS elevations are typically reported as departures from the long-term ELA₀. The ELA₀ for Brewster Glacier was initially estimated using AAR, assuming the ratio of the accumulation area to the full glacier area was equal to 2:1 or 0.66. As the EOSS programme matured, the initial ELA₀ was adjusted by regressing annual departures from ELA₀ on Brewster Glacier against the annual mean departure from ELA₀ measured across the Southern Alps (ELA_Alps). The constant of the regression is used to adjust ELA₀ assuming a zero average mass balance for the Southern Alps (Chinn and others, 2012). This adjustment has been applied on two occasions to Brewster Glacier, resulting in ELA₀ changing from 1880 to 1905 m a.s.l. in 1999 (AAR₀ = 0.6) and to ELA₀ = 1935 m a.s.l. (AAR₀ = 0.5) in 2003. As noted by Willsman and others (2015, p. 16) estimates of ELA₀ are ‘dependent on the accuracy of the digitised accumulation and/or ablation areas, glacier area and its area elevation curve’. In other words, the glacier outline and DEM being used control the relationship between ELA₀ and AAR₀. For example, applying the 1:1 ratio (AAR = 0.5) using the dataset in this study results in ELA₀ = 1918 m a.s.l., 17 m below the ELA₀ reported for EOSS computations. The calculation of AAR is also subject to variability depending on the approach used to determine its value. Two-dimensional mapping of the equilibrium line in b _eq yields AAR_eq = 0.4 for ELA_eq = 1923 m a.s.l., while a 1-D approach using an area elevation curve yields an AAR of 0.46, which is closer to the EOSS-determined AAR₀. Clearly, these examples demonstrate the obvious, which is the calculation of ELA and/or AAR is sensitive to the method chosen to do it (e.g. Braithwaite, 2015).

6.5. Theoretical relationship between ELA, AAR and mass balance

The relationship between ELA, AAR and mass balance has been central to revealing linkages between glacier behaviour and climate, and has been used effectively in efforts to monitor changes in mass balance using remote-sensing methods. To determine mass balance from ELA and/or AAR using remote sensing it is important that the effects of glacier hypsometry are taken into account. The area-altitude distribution of a glacier can distort the relationship between mass balance and ELA or AAR, and may result in the commonly-used linear relationship being valid only within certain elevation ranges (Khalsa and others, 2004). Our glacier in equilibrium $(b_{a_{eq}})$ map (Fig. 12) provides us with an opportunity to explore and characterize the shape of such relations for Brewster Glacier. Figure 13 shows that despite inflexion points the AAR relationship is close to linear over the annual mass-balance range (−3000 to 2000 mm w.e. When AAR exceeds 0.95 (glacier nearly fully snow covered at the end of the ablation season), it is clear that the linear relationship would underestimate the very positive balance. The relationship between mass balance and ELA is also linear until it reaches an altitude of ~2050 m a.s.l. At and above this altitude, ELA would be located on the upper, steep slopes and if it were to continue to increase in elevation the change in mass balance would be small. All data points fit the theoretical shapes well, except in 2011 when the mass balance is similar to 2008 but ELA is much higher. It is possible that in this very negative year, the regime of ablation departed from its typical shape predicted by MSBI _CK using linear regression. This is supported by the relatively low correlation between stake ablation and MSBI _CK for this year (R ² = 0.62). Nonetheless, the relationships shown in Figure 13 are important for our understanding of the behaviour of Brewster Glacier and are likely to be useful in reconstructing mass balance using remote-sensing methods. They also shed new light on whether it is always appropriate to assume linearity for such relationships.

Fig. 13. Relationships between mass balance and (a) AAR and (b) ELA for Brewster Glacier. The solid line is derived using $b_{a_{eq}}$ . The envelope corresponds to uncertainties associated with B _{a t} . The dots indicate the values derived from the modelling of annual balance for each year, with error bars showing $\sigma _{ELA_t}$ . The dotted lines in (b) show the expected magnitude of $\sigma _{ELA_t}$ .

6.6. Mass-balance gradient

The mass-balance gradient, defined as the rate of change of mass balance with altitude, is an important parameter in determining mass balance using remote-sensing methods (e.g. Rabatel and others, 2005; Chinn and others, 2012). Extrapolation of a known mass-balance gradient to neighbouring regions characterised by similar climatic conditions enables mass balance to be determined in areas where no glaciological field data are available. This approach is used by Chinn and others (2012) to estimate ice volume changes in the Southern Alps using the EOSS record, with mass-balance gradients determined from glaciological measurements on three glaciers, including Brewster Glacier. Mass-balance gradients in the ablation area of glaciers in the Southern Alps are assumed to be twice those of the accumulation area by Chinn and others (2012), while higher values are assigned to the ‘wet’ western glaciers, and lower values to the ‘dry’ eastern glaciers. Following this approach, mass gradients west (east) of the Main Divide of the Southern Alps are defined as 25.8 (15) and 12.9 (7.5) mm w.e. m⁻¹ for the ablation and accumulation areas, respectively. For Brewster Glacier, the average mass-balance gradients are reported as 19.9 and 6.8 mm w.e. m⁻¹ for the ablation and accumulation areas, respectively.

In this context, it is worth examining the typical altitudinal mass-balance gradient from $b_{a_{eq}}$ . Figure 14 shows the nature of the relationship between mass balance and elevation, as well as the derived gradient. Although the values are in general agreement with the estimates derived by Stumm (2011), Figure 14 reveals some additional information. The mass-balance gradient peaks twice, namely in the lower ablation zone and around ELA_eq, with the gradient approaching or exceeding 25 mm w.e. m⁻¹. The elevation band between 1840 and 1920 m a.s.l. has a lower gradient, which is similar to the gradient of 5 mm w.e. m⁻¹ above 2000 m a.s.l. On average, the altitudinal gradients of $b_{a_{eq}}$ are equal to 14.5 and 7.4 mm w.e. m⁻¹ in the ablation and accumulation zones, respectively (i.e. upward and downward of ELA_eq, respectively). These depart slightly from the values reported by Chinn and others (2012) for Brewster Glacier but differ from those expected for a glacier that is located west of the Main Divide. In fact, they are quite similar to the average mass-balance gradients estimated for Tasman Glacier and the mean values used for the ‘dry’ eastern index glaciers. This raises the question as to whether it is appropriate to simply partition glaciers into two mass-balance gradient categories based on whether they are west or east of the Main Divide (e.g. Chinn and others, 2012). There is also considerable variability in the average annual mass-balance gradient on Brewster Glacier (Table 3), as well as at the ELA and in the ablation/accumulation zones. The largest variation in the mass-balance gradient is in the ablation zone, which is consistent with the finding that summer balance tends to have the largest influence on annual mass balance. Beyond documenting this record, the observed variability is likely to be useful in determining uncertainties in deriving mass balance using methods such as proposed by Rabatel and others (2005).

Fig. 14. Altitudinal mass balance and derived gradient of Brewster Glacier. The envelope corresponds to uncertainties associated with B _{a t} . The dotted lines depict the spread of $b_{a_{eq}}$ along altitudinal contours.

Table 3. Mass-balance gradients for each year

The values in bold show the mean and standard deviation of the mass-balance gradients.

6.7. Limitations of the method

Our proposed interpolation scheme makes the assumption that indices of snow depth and ablation do not change over time. This is justified in this study due to the relatively short period of analysis, but the assumption would need to be scrutinised over a longer time period. For example, changes in the surface elevation and topography of Brewster Glacier over time or a shift in climate are likely to lead to changes in the regimes of accumulation and/or ablation. In some respects, our approach is similar to that proposed by Lliboutry (1974), who developed a linear model to calculate mass balance using site (spatial term) and year (temporal term) as independent variables. By using all available data over 16 consecutive years (1956–72) from the ablation area of Glacier de Saint-Sorlin (French Alps) to regress both terms, Lliboutry (1974) mitigated problems associated with observational error from any particular year. The multi-variate statistical method has also been used in conjunction with geodetic measurements to enable estimates of mass balance (Thibert and Vincent, 2009). Eckert and others (2011) further refined the method by relaxing the assumption of stationarity, enabling temporal trends to be characterised in a mass-balance series obtained from Glacier de Sarennes (French Alps). However, this approach does not explicitly account for spatial variability in the mass-balance record. By focusing on both the spatial and temporal structure of mass-balance variability using a geostatistical method, while utilising the entire dataset available, our method to some extent fills the gap between the purely statistical approaches of Lliboutry (1974) and Eckert and others (2011) and the more traditional semi-empirical methods of interpolating (extrapolating) glaciological observations using altitude. We do this by using our indices to represent to some degree the spatial term proposed by Lliboutry (1974), while the temporal term is captured by the slope and offset of the linear regression between our glaciological measurements and the indices. Further observations spread over a longer timeframe will enable the indices to be refined spatially, as well as to reduce their standard error under the assumption of stationarity. In the future, new indices could be considered for different timeframes to explore temporal variations in the spatial structure of accumulation and ablation. Improvements to our mass- balance estimates could also be achieved through validation using geodetic mass-balance estimates (e.g. Zemp and others, 2013) and/or homogenization of the record using a distributed energy-balance model (e.g. Huss and others, 2009). The latter would overcome the deficiencies associated with mass balance being evaluated between the dates of successive field surveys, rather than an evaluation between fixed dates (e.g. a hydrological year).

7. CONCLUSION

We have reported on the glaciological measurements used to reconstruct mass balance on Brewster Glacier. The efforts to gather glaciological observations on Brewster Glacier have been considerable, with a large number of people contributing in different ways since 2004. Despite the logistical challenges, it is the longest in situ record of a glacier's mass change ever obtained in the Southern Alps and as a result, has become an important benchmark for a range of glaciological and meteorological studies. To reduce spatial and temporal artefacts introduced by changes to the measurement regime as the mass-balance programme evolved, we re-analysed all the observational data and developed a new geostatistical method to calculate mass balance. We have shown that the spatial structure of winter and summer balance over Brewster Glacier for any individual mass-balance year contains sufficient similarity to other years to enable indices to be determined for both. The newly derived indices are more powerful than elevation in revealing the spatial structure of winter and summer balance. Changes in winter and summer mass balance were obtained directly by calibrating the indices with data points available for any given year, which helped distribute errors associated with potential artefacts introduced from outliers from any individual mass-balance year. The method also enabled a predictor for annual balance to be resolved, which was used specifically for years where values of summer balance were not available.

The average glacier-wide surface winter, summer and annual mass balances for the period 2005–15 were 2484, −2586 and −102 mm w.e., respectively. The variability in summer balance was greater than in the winter balance, with the strong relationship between summer and annual balance (R ² = 0.87, P-value < 0.01) suggesting ablation processes exert the largest control on annual balance. On the whole, there is good agreement between our ELA and AAR values and those derived from the EOSS programme (Chinn and others, 2005, 2012), though discrepancies in some years could not be fully accounted for. A mass-balance map of Brewster Glacier in an equilibrium state, which by definition has a glacier-wide mass balance equal to zero, was used to calculate values of ELA (1923 ± 25 m a.s.l.) and AAR (0.40) representative of the observational period. The relationships between mass balance and ELA/AAR were explored, demonstrating they are mostly linear, which provides a foundation for reconstructing mass balance using remote-sensing methods. On average, the mass-balance gradients were found to be equal to 14.5 and 7.4 mm w.e. m⁻¹ in the ablation and accumulation zones, respectively, but there is considerable variability from year to year, especially in the ablation zone. We also observed variability in the mass-balance gradients between different elevation bands, clearly demonstrating that care must be taken when prescribing fixed mass-balance gradients to estimate glacier volume changes in the Southern Alps.

The geostatistical approach used in this research provides a robust and transparent method to assess the spatial and temporal variability of mass balance, which should allow the data from Brewster Glacier to be easily integrated into global glacier monitoring databases. This mass-balance record provides a baseline for future research, and has already been useful in efforts to build new linkages between large-scale atmospheric circulation and glacier behaviour in the Southern Alps, as well as to reconstruct glacier behaviour over the EOSS period (1977-present) using a remote-sensing method (e.g. Sirguey and others, 2016). Because the method takes advantage of the entire glaciological record to help refine the predictors of mass balance (defined as indices), it will continue to be improved over time, which will allow issues of non-stationarity to also be assessed. In the short term, observational data from the upper accumulation area will be gathered to help reduce uncertainties in the calculation of mass balance. A geodetic survey of mass balance and/or homogenization of the observed mass-balance record using a distributed energy-balance model are the next logical steps to further scrutinise the newly constructed record.