Reconstructing the mass balance of Brewster Glacier , New Zealand , using MODIS-derived glacier-wide albedo

In New Zealand, direct measurements of mass balance are sparse due to the inaccessibility of glaciers in the Southern Alps and the logistical difficulties associated with maintaining a mass balance record. In order to explore the benefit of remotely sensed imaging to monitor mass balance in the Southern Alps, this research assesses the relationship between measurements of glacier surface albedo derived from Moderate Resolution Imaging Spectroradiometer (MODIS) and mass balance observations using the glaciological method on Brewster Glacier over the 2005–2013 period. We confirm that minimum glacier-wide albedo is a reliable predictor for annual mass balance in this maritime environment (R2 = 0.93). Furthermore, we show that regular monitoring of glacier-wide albedo enables a new metric of winter accumulation to be derived, namely the cumulative winter albedo, which is found to correlate strongly with winter mass balance (R2 = 0.88), thus enabling the reconstruction of separate winter and summer mass balance records. This allows the mass balance record for Brewster Glacier to be extended back to the start of MODIS observations in 2000 and to confirm that the annual balance of Brewster Glacier is largely controlled by summer balance (R2= 92 %). An application of the extended record is proposed whereby the relationship between mass balance and the photographic record of the end-of-summer snowline altitude is assessed. This allowed the annual balance record of Brewster Glacier to be reconstructed over the period 1977–2013, thus providing the longest record of mass balance for a glacier in New Zealand. Over the 37-year period, our results show that Brewster Glacier gained a significant mass of up to 14.5± 2.7 m w.e. by 2007. This gain was offset by a marked shift toward negative balances after 2008, yielding a loss of 5.1± 1.2 m w.e., or 35 % of the gain accumulated over the previous 30 years. The good correspondence between mass balance of Brewster Glacier and the phase of the Pacific (Inter-)Decadal Oscillation (PDO/IPO), associated with the fast terminus retreat observed between 1978 and 1998, strongly suggests that the observed mass gain of Brewster Glacier since 1977 is only offsetting a longer sequence of dominantly negative balances.

Abstract.In New Zealand, direct measurements of mass balance are sparse due to the inaccessibility of glaciers in the Southern Alps and the logistical difficulties associated with maintaining a mass balance record.In order to explore the benefit of remotely sensed imaging to monitor mass balance in the Southern Alps, this research assesses the relationship between measurements of glacier surface albedo derived from Moderate Resolution Imaging Spectroradiometer (MODIS) and mass balance observations using the glaciological method on Brewster Glacier over the 2005-2013 period.We confirm that minimum glacier-wide albedo is a reliable predictor for annual mass balance in this maritime environment (R 2 = 0.93).Furthermore, we show that regular monitoring of glacier-wide albedo enables a new metric of winter accumulation to be derived, namely the cumulative winter albedo, which is found to correlate strongly with winter mass balance (R 2 = 0.88), thus enabling the reconstruction of separate winter and summer mass balance records.This allows the mass balance record for Brewster Glacier to be extended back to the start of MODIS observations in 2000 and to confirm that the annual balance of Brewster Glacier is largely controlled by summer balance (R 2 = 92 %).An application of the extended record is proposed whereby the relationship between mass balance and the photographic record of the end-of-summer snowline altitude is assessed.This allowed the annual balance record of Brewster Glacier to be reconstructed over the period 1977-2013, thus providing the longest record of mass balance for a glacier in New Zealand.
Over the 37-year period, our results show that Brewster Glacier gained a significant mass of up to 14.5 ± 2.7 m w.e. by 2007.This gain was offset by a marked shift toward negative balances after 2008, yielding a loss of 5.1 ± 1.2 m w.e., or 35 % of the gain accumulated over the previous 30 years.The good correspondence between mass balance of Brewster Glacier and the phase of the Pacific (Inter-)Decadal Oscillation (PDO/IPO), associated with the fast terminus retreat observed between 1978 and 1998, strongly suggests that the observed mass gain of Brewster Glacier since 1977 is only offsetting a longer sequence of dominantly negative balances.

Introduction
Many of the world's small mountain glaciers are losing mass at an accelerated rate in response to climate warming over the last century (Zemp et al., 2015).Annual fluctuations in glacier volume are best quantified using the well-established, but labour-intensive glaciological method.These measurements represent the direct response of the glacier to meteorological conditions that control ablation and accumulation, with long-term mass balance trends providing a valuable climate signal (Six and Vincent, 2014).At a broader scale, mass balance records from multiple glaciers are useful for estimating the contribution of glacial melt to sea level rise (Kaser et al., 2006;Bahr et al., 2009).P. Sirguey et al.: Reconstructing the mass balance of Brewster glacier Despite the scientific importance of glacier mass balance records, the global distribution of glaciers with direct mass balance monitoring programmes is strongly biased towards the Northern Hemisphere (Zemp et al., 2015).The glaciers within the Southern Alps of New Zealand are among the few temperate glaciers in the Southern Hemisphere and are, therefore, regarded as a single glacier region with global significance (Radić and Hock, 2010;Pfeffer et al., 2014).However, only one New Zealand glacier, Brewster Glacier, has a current and continuous direct mass balance record exceeding 10 years (2005 to present) (Cullen et al., 2016).In order to extend the temporal window of this record, this paper evaluates an alternative method to estimate mass balance by monitoring the evolution of glacier surface albedo derived from the Moderate Resolution Imaging Spectroradiometer (MODIS).Given the logistical difficulties associated with initiating and maintaining mass balance monitoring programmes in alpine environments, such a method based on remotely sensed imagery may prove useful in helping to reconstruct mass balance records for other New Zealand glaciers.
In lieu of direct glaciological monitoring, the equilibrium line altitude (ELA) is regarded as a good proxy for annual balance (Braithwaite, 1984;Rabatel et al., 2005).Accordingly, ELA measurements from remotely sensed data have been used as a substitute for direct mass balance monitoring (Chinn et al., 2005).This is the approach currently used in New Zealand to infer the mass balance variability of 50 "index" glaciers situated in the Southern Alps, including Brewster Glacier (Chinn et al., 2012).For these glaciers, the position of the end-of-summer snowline (EOSS, an approximation of the ELA) is identified from oblique photographs obtained from annual aerial surveys conducted since 1977.Positive or negative mass balances are inferred by comparing the annual EOSS altitude (hereafter referred to as SLA) to the steady-state ELA 0 with assumptions on regional mass balance gradients at the ELA (Rabatel et al., 2008;Chinn et al., 2012).While the EOSS monitoring programme has provided an exceptional amount of data on the behaviour and climate response of New Zealand's glaciers (e.g.Fitzharris et al., 1997;Clare et al., 2002), the reconstruction of mass balance for a single glacier from this record has not yet been attempted.
In relation to potential mass balance proxies, previous studies have reported on linkages between satellite-derived measurements of surface albedo and annual balance (Greuell et al., 2007;Dumont et al., 2012;Brun et al., 2015).In particular, Dumont et al. (2012) reported on a strong correlation between annual balance and annual glacier-wide albedo minima retrieved from MODIS data for St Sorlin Glacier in the French Alps.Although albedo minima are not regarded as a "conventional" proxy for mass balance, the annual minimum glacier-wide albedo (reached at the end of the ablation season) is of significance because it measures the relative share of exposed ice in the ablation area and snow and firn in the accumulation area.It is thus closely related to the accumulation area ratio (AAR), itself being related to annual balance (Racoviteanu et al., 2008;Dyurgerov et al., 2009;Bahr et al., 2009).The validity of this approach to estimating mass balance was further assessed by Brun et al. (2015) on two Himalayan glaciers with mixed results.The good relationship between minimum glacier-wide albedo and observed annual balance for Chhota Shigri Glacier was used to reconstruct a mass balance record over the period of MODIS observations.The less reliable relationship for Mera Glacier was attributed to a less marked and poorly resolved albedo cycle due to the persistent cloud cover during the Indian monsoon.These mixed results justify further investigation into the suitability of this method for use in the New Zealand context to extend the mass balance record of Brewster Glacier beyond glaciological observations, as well as assess whether monitoring albedo on New Zealand glaciers can provide insights about glacier fluctuations.
In this context, this study first assesses the accuracy of MODIS-derived albedo measurements against field albedo data.Based on the method originally proposed by Dumont et al. (2012) we assess the relationships between albedo measurements derived from MODIS data and the annual balance of Brewster Glacier.This allows the mass balance series to be extended over the time period 2000-2013 when MODIS/TERRA images are available.We further develop the method by showing how the evolution of glacier-wide albedo over a mass balance year can be used to reconstruct both winter and summer mass balance records.This provides 5 additional years (2000)(2001)(2002)(2003)(2004) of seasonal mass balance data that may provide valuable insights into the accumulation and ablation processes governing the evolution of Brewster Glacier (Cullen et al., 2016).Finally, an application of this method is presented whereby the extended mass balance series is used to reconstruct the mass balance record through the full length of the EOSS programme .
The next section provides a description of Brewster Glacier, the observed mass balance record, the field albedo measurements, and the required MODIS data products.In the third section, the albedo retrieval algorithm and the mass balance reconstruction method are described.The fourth section considers the results of the MODIS accuracy assessment and evaluates the relationships between MODIS-derived albedo, mass balance, and SLA.The reconstructed mass balance records are also presented.The final section discusses the reconstructed sequence of mass balance of Brewster Glacier and its links to regional climate in the Pacific region.
2 Site description and data sets

Brewster Glacier
Brewster Glacier is a small maritime glacier, predominately free of debris, with an area of ca. 2 km 2 located on the western side of the Southern Alps in Mount Aspiring National Park, New Zealand (44.08 • S, 169.43 • E) (Fig. 1).The glacier descends a distance of ca.2.5 km from a maximum elevation of 2400 m a.s.l. to an elevation of 1700 m a.s.l.At higher altitudes (> 2000 m a.s.l.), the glacier has a southwesterly aspect and a mean slope of 31 • while the lower portion of the glacier has a southerly aspect and a mean slope of 10 • .Brewster Glacier is one of the most well-researched New Zealand glaciers, with recent studies focusing on the atmospheric/meteorological controls on the surface energy and mass balance (Gillett and Cullen, 2011;Conway et al., 2015;Conway andCullen, 2013, 2016;Cullen and Conway, 2015), and the climate sensitivity of Brewster Glacier (Anderson et al., 2010).Brewster Glacier belongs to the 50 index glaciers monitored by the EOSS survey programme with a record of the end of summer snowline altitude (SLA) since 1978 (Chinn et al., 2005;Willsman et al., 2015).The current in situ mass balance monitoring programme began in February 2004, making it the longest mass balance record using the glaciological method for a New Zealand glacier.This mass balance record has recently been reanalysed and consolidated by Cullen et al. (2016).The overlap between glaciological mass balance, MODIS, and EOSS records provides an opportunity to reconstruct the mass balance record of Brewster Glacier back to 1977 when the EOSS programme started and thus improve knowledge about its long-term evolution.

Glacier outline and surface topography
For the purpose of this study, a 15 m-resolution DEM interpolated from 20 m contours derived from a photogrammetric aerial survey in 1986 (NZSoSDEM, Columbus et al., 2011) was used for all calculations involving height due to its consistent representation of the surface.The glacier boundary 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 referencesurface mass balance (Elsberg et al., 2001).

Glaciological observations
The mass balance programme on Brewster Glacier was initiated in February 2004 whereby the traditional glaciological method is used to measure accumulation directly using snow pits and probing, whereas stakes are used to measure surface ablation.This record of direct glaciological observations was reanalysed using a geostatistical approach described by Cullen et al. (2016), yielding a consolidated and homogenized time series of direct mass balance observations (Fig. 2).Presently, the winter balance is determined at the end of the winter (early to mid-November) from geolocated snow depth measurement (snow probing) and snow density measurements from at least one snow pit dug on the central flow line.Snow and ice ablation is determined by measuring changes in the height of a network of PVC stakes drilled into the glacier using a Heucke steam ice drill (Heucke, 1999) during the accumulation survey campaign, thus allowing summer balance to be derived and, in turn, the annual balance.

Ground albedo measurements
Albedo measurements retrieved from MODIS images were first validated by comparison to a 22-month albedo record from an automatic weather station (AWS) installed in the ablation zone of Brewster Glacier from 25 October 2010 to  (Cullen and Conway, 2015).Note that an upper threshold of 0.95 was applied when SW↓ was compromised by snow covering the upper radiation sensors during heavy snowfall (35 occurrences out of 678 days).
1 September 2012 (Cullen and Conway, 2015).The AWS was installed on the central flow line of the glacier at an elevation of 1760 m a.s.l.At this location, the glacier surface has a south-south-westerly aspect and a mean slope of 6 • .Incoming and outgoing shortwave radiation (SW↓ and SW↑ respectively) were measured every 30 s with a Kipp and Zonen CNR4 net radiometer and stored as 30 min averages.The net radiometer was mounted at an initial height of 2.1 m above the surface and changes in the surface height of the glacier were accounted for by periodic raising or lowering of the AWS (mounted on an aluminium pole drilled into the glacier) as required.At a height of 2.1 m, the lower sensor/downward-facing pyranometer captured a footprint of ca.200 m 2 .Between January and March 2011, high ablation rates caused tilting of the AWS mast which can yield errors in measured albedo due to the underestimation of SW↓ (Bogren et al., 2016;Wang et al., 2016).To account for this as well as the effect of surface slope, incoming shortwave radiation data were adjusted using the tilt and slope correction method developed by van As (2011).Albedo data were obtained from the ratio between SW↑ and SW↓, and processed into hourly accumulated albedo, that is, the ratio of accumulated SW↑ and SW↓ over an hour (van den Broeke et al., 2004), for comparison with MODIS-derived albedo measurements.Although there is no standard that addresses the net measurement accuracy of pyranometers, CNR4 are First Class ISO9060 pyranometers including calibration uncertainty (< 1 %), worst case directional response (< 2 %), temperature response (< 4 %), annual non-stability (< 1 %), and tilt error (< 1 %) (Kipp and Zonen, 2014).This led to a net uncertainty of incoming and outgoing solar radiation measurements estimated at ±5 %, yielding in turn an uncertainty of ±7 % on the resulting albedo value.The 22-month record is illustrated in Fig. 3 by the time series of daily accumulated albedo as defined by van der Broeke et al. (2004).

EOSS programme
Although records of SLA on selected New Zealand glaciers started in 1977, the first aerial photo of Brewster Glacier was captured in 1978 (Willsman et al., 2015).The long-term equilibrium line (ELA 0 ) of Brewster Glacier from which the departure of the SLA for each year is calculated was estab-lished at 1935 m a.s.l.(Willsman et al., 2015).(Chinn et al., 2005).The linear relationship for Brewster (see Willsman et al., 2015) enables missing dates to be estimated from the SLA Alps record of the remaining index glaciers.For the years 1979, 1990, and 1991, when no flight was conducted and for which ground records only existed for five, two, and one glaciers respectively, SLA Alps was first determined using its relationship with the SLA of each glacier on record that year.This was done by finding the weighted least-squares solution of the system of linear equations formed by the relationship associated with each glacier, with the corresponding coefficients of determination being used as weights.Figure 4 illustrates the observed and estimated SLA departures for Brewster Glacier.The uncertainty associated with the SLA record is difficult to characterize because it applies to various methodologies, each affected by interpretation and subjectivity (Willsman et al., 2015).For example, the SLA is often determined by applying an interpolation method, which involves arranging all photographs into a sequence based on increasing area of snow cover (descending order of SLAs)."The photograph being interpolated is then carefully compared and inserted into its appropriate place in the sequence" (Willsman et al., 2015, p. 15).Despite such subjectivity, Chinn et al. (2012) estimated the uncertainty on the mean annual SLA departure from the 50 index glacier to be ±5 m at the 95 % confidence level.This would correspond to an estimated standard error for an individual SLA to be ±18 m.In the context of this study, we used an uncertainty of ±30 m, which we believe to be conservative for this data set.

MODIS products
The success of the mass balance reconstruction is dependent on our ability to derive accurate measurements of glacier surface albedo from MODIS data.As Brewster Glacier is a  (Willsman et al., 2015).
small mountain glacier surrounded by steep, rugged topography, the retrieval of albedo from satellite data is complicated by the relatively small size of the glacier combined with the effects of a steep and highly reflective topography on the radiative budget.Commonly used MODIS albedo products such as MOD10 daily snow product or MOD43 16day albedo product are not appropriate for Brewster Glacier because the spatial resolution (500 m) is too coarse to resolve albedo variability across the glacier surface.Furthermore, Brun et al. (2015) show that the alternative MODIm-Lab albedo product chosen for this study (see Sect. 3.1) significantly outperformed the MOD10 albedo product.
The MODImLab albedo retrieval algorithm uses MODIS Level-1B Swath products (MOD02QKM, MOD02HKM, MOD021KM, MOD03) that provide radiance measurements at the top of atmosphere (TOA) in the 36 reflective and emissive spectral bands of the MODIS sensor, as well as relative geometry parameters such as solar and sensor zenith (θ s , θ v respectively) and solar and sensor azimuth.For this study, all MODIS/TERRA Level 1B C5 products acquired between 24 February 2000 and 31 December 2013 over Brewster Glacier were obtained.This provided 6871 granules captured on 4961 unique days, which reduced to 6187 images after mosaicking subsequent granules.842 images were available over the 22-month AWS observation period for the purpose of validating the MODIS-derived albedo.

MODImLab snow and ice albedo product
MODImLab is a MATLAB toolbox allowing MODIS L1B swath data products to be processed to map subpixel snow fraction (Sirguey et al., 2009;Charrois et al., 2013), snow and ice albedo (Dumont et al., 2011(Dumont et al., , 2012)), and snow specific surface area (Mary et al., 2013).It implements an image fusion approach that improves mapping in steep terrain by combining MOD02QKM bands 1 and 2 (250 m resolution) with MOD02HKM bands 3 to 7 (500 m resolution) to produce seven spectral bands at a spatial resolution of 250 m (Sirguey et al., 2008).At its core is the ATOPCOR atmospheric and topographic correction module that corrects images of TOA radiance into maps of ground surface reflectance (Sirguey et al., 2009).It accounts for multiple reflections of radiation between the glacier surface and surrounding topography (Sirguey, 2009), which has been shown to improve mapping of snow and albedo in mountainous terrain (Dumont et al., 2011).White-sky (WS) and bluesky (BS) snow and ice albedo are retrieved using a narrowband to broadband conversion based on a snow radiative transfer model and can benefit from the implementation of a snow anisotropic reflectance model (Dumont et al., 2010(Dumont et al., , 2011)).MODImLab also implements a cloud detection algorithm based on a selection of tests largely inspired by the MOD35 product (Ackerman et al., 2006) and applied to the MODIS reflective and emissive bands from the MOD021KM swath L1B products.Dumont et al. (2012) estimated the accuracy of MODImLab albedo output (α MOD ) to be better than ±10 % under clear-sky conditions.As recommended by Dumont et al. (2012), we used WS albedo due to its lower sensitivity to seasonal variations in illumination.Figure 5 illustrates an output of α MOD retrieved with the anisotropic reflectance model.

Characterization of α MOD accuracy
Of the 842 images acquired during the ground albedo measurements campaign, 497 were excluded because clouds obscured the AWS pixel according to the MODImLab cloud mask (n = 452) or the measured albedo was considered to be invalid (α AWS > 1, n = 45).The remaining 345 images were used to characterize the performance of the albedo retrieval algorithm.MODIS-derived albedo measurements α MOD were extracted from the pixel corresponding to the AWS site, and compared to the hourly accumulated albedo measurements α AWS determined from both 30 min observations that bracket the satellite overpass time.In order to facilitate the production of a usable 14-year time series of α MOD , we assessed the agreement between α MOD and α AWS under different scenarios of cloud, illumination/viewing geometry, and using the isotropic or anisotropic reflectance model (Dumont et al., 2011).This yielded criteria and thresholds for the automatic exclusion of images with poor-quality albedo as a first filtering step (see Sect. 3.3.1).
For example, large errors in α MOD can be caused by nearby/misclassified clouds in the MODImLab cloud mask that compromise the clear-sky hypothesis of ATOP-COR (Brun et al., 2015).Although systematic visual inspection of processed imagery can be used to exclude obvious un- reliable data points from the time series, we aimed at avoiding the need for manual interventions.We, therefore, characterized the effect of relying exclusively on the MODImLab cloud mask by examining the AWS SW↓ record to identify a conservative subset of the 345 images for which clear-sky conditions were guaranteed at the time of MODIS observations.The subset of 263 images was determined by visually assessing the shape of the diurnal cycle of measured SW↓ compared to that predicted by integrating the SPCTRAL2 clear-sky solar spectral model of Bird and Riordan (1986) over the range of wavelengths captured by the pyranometer (i.e.300-2800 nm). Figure 6 illustrates this selection by displaying successive days of measured SW↓ compared to modelled clear-sky solar irradiance in both summer and winter time.
Complex illumination/viewing geometry may also compromise the quality of the α MOD maps.Wolfe et al. (1998) document that MODIS pixel size increases 2-to 5-fold from θ v = 45 • to a maximum of 66 • .The pronounced panoramic effect of MODIS scans thus potentially yields increasing spectral mixtures from larger pixel footprints away from nadir that can compromise the albedo retrieval algorithm.Large solar and sensor zenith values may also exacerbate the effect of snow/ice anisotropic reflectance, although Dumont et al. (2011Dumont et al. ( , 2012) ) showed that the use of the anisotropic model can mitigate this issue and improve the albedo retrieval.The sensitivity of the albedo retrieval to θ v and θ s , and the effectiveness of the anisotropic reflectance correction of snow and ice was assessed by examining the error α = α MOD − α AWS as a function of the illumination/viewing geometry.Several performance metrics were used to characterize α , namely the root mean square deviation (RMSD), the mean bias, the mean absolute error (MAE), the Pearson's coefficient of determination (R 2 ), and linear regressions.

MODIS-derived mass balance
Dumont et al. ( 2012) and Brun et al. (2015) give evidence that the evolution of albedo at the surface of a glacier contains information about mass balance.This can be explained first by the fact that albedo often controls a significant share of the energy budget of glaciers.Secondly, if considering the glacier surface to be split into two contrasted targets, namely highly reflective snow in the accumulation area ( α ≈ 0.8), and less reflective exposed ice (α ≈ 0.4) in the ablation area, then glacier-wide albedo is a case of linear spectral mixing and can be shown to approach the AAR, itself often a predictor of mass balance (Dyurgerov et al., 2009).In this context, the albedo averaged across the surface of a glacier, hereby named glacier-wide albedo α, reaches a minimum during the ablation season of year yr noted α min yr that can be linearly correlated to annual mass balance B a yr .Taking advantage of the high temporal resolution of MODIS observations Dumont et al. (2012) and Brun et al. (2015) resolved α min yr over several years and built a linear model relating to the observed mass balance of several glaciers.These regression models are in turn used to estimate mass balance over the length of the MODIS record.In this study we test this method on Brewster Glacier, but also further explore the time series of glacierwide albedo to find predictors of seasonal (i.e.summer and winter) mass balance.

Time series of MODIS-derived glacier-wide albedo
After assessing the accuracy of α MOD from Brewster Glacier (see Sect. 4.1), we applied the albedo retrieval algorithm with anisotropic correction to all MODIS Level 1B products available between February 2000 and December 2013.Only maps with no clouds within the glacier margin according to MODImLab were retained.Although the best accuracy of albedo was found when θ v < 30 • , images acquired up to θ v < 45 • were included as this nearly doubled the number of available albedo images.Despite the lower performance achieved with θ v < 45 • , this was believed to be an acceptable compromise to better resolve the temporal variability of albedo, in particular when it reached its minimum.The complete albedo record (2000-2013) consisted of 913 maps of albedo.The glacier-wide albedo α(t) was calculated for each map acquired at time t from the 26 pixels shown in Fig. 5.In order to mitigate remaining noise, measurement errors, and high frequency variability in the albedo signal, the time series was smoothed with a 3 period rolling average.Because of the unevenly spaced time series, this method ensured that all smoothed values benefited from a similar statistical averaging as opposed to a fixed time rolling average, which would have involved a variable number of samples for each smoothed value.

Cumulative winter albedo
Resolving winter and summer mass balance on any given year yr (B w yr and B s yr respectively) helps to inform us about processes controlling the variability in annual balance on Brewster Glacier (Cullen et al., 2016).Therefore, giving consideration to the evolution of α(t) over a mass balance year can potentially inform us about the variability in winter snow accumulation.In this context, we hypothesize that winter accumulation is positively correlated to the time during which the glacier surface is fully covered by fresh snow or snow in the early stages of metamorphism.In addition to contributing directly to the mass gain, fresh snow highly reflects shortwave radiation and, therefore, also minimizes ablation.Thus, long periods of highly reflective full snow cover capture several situations contributing positively to the annual balance (e.g.large snowfall, frequent snow events, long accumulation season, and delayed exposition of ice).We define the cumulative winter albedo A w yr of accumulation season yr as the integral of α(t) between 1 March and 31 October of year yr (to match the measurement period of B w ) when α(t) exceeds the 72 % threshold.We found that this threshold yielded the highest correlation with B w and is taken to be representative of fresh or recent snow fully covering the glacier surface.For the purpose of computing A w yr , α(t) was linearly interpolated for all daily time steps.The integral has the benefit of more heavily weighting the days with high albedo so that the days exhibiting a more reflective surface, indicative of new accumulation and/or reduced ablation, contribute more to A w yr .

MODIS-derived mass balance models
We investigated linear regression models between values derived from α(t) (i.e.α min yr and A w yr ) and the annual (B a yr ) and seasonal (B w yr , B s yr ) mass balances of Brewster Glacier estimated from glaciological observations (Cullen et al., 2016).Two models were considered.
1. Model M1 (similar to Dumont et al., 2012 andBrun et al., 2015): 2. Model M2: where s i and o i are the slope and offset of each linear regression respectively.
The standard error for the mass balance data has been characterized by Cullen et al. (2016).The uncertainty associated with α min yr is obtained as the standard error RMSD √ 26 , where RMSD is obtained from the accuracy assessment of α MOD .The uncertainty of A W yr is related to the number of days N during which α (t) > 72 % and is given by RMSD √ N 26 .Since there is uncertainty in all variables, rigorous linear regression models were calculated using the unified regression equations of York et al. (2004) instead of ordinary least squares.The standard error in the regression parameters σ slope could then be used to derive the standard error of prediction for new observations (Montgomery and Runger, 2011;421-423).
The better performance of Model M1 at estimating annual balance (see Sect. 3.3) justified its use to extend the mass balance record over the period of MODIS observations.The coefficient of determination R 2 and standard error of prediction also indicated a better estimation of B w compared to B s in Model M2.Therefore, in order to provide estimates

Mass balance reconstruction
In addition to assessing whether the mass balance of a New Zealand glacier can be estimated from monitoring its surface albedo, the proposed method permits the mass balance record of Brewster Glacier to be extended beyond the glaciological record with MODIS-derived estimates of annual balance B MODIS a yr , as well as seasonal balances for the 2000-2004 period.These five additional estimates from MODIS were used to supplement the glaciological estimates of the 2005-2013 period yielding a 50 % increase of the record length.As an application of this extended record, we assessed the relationship between the departure of SLA from ELA 0 determined by the EOSS programme and the extended mass balance record in order to reconstruct the annual balance over the length of the EOSS programme (1977-present).
In this context, Cullen et al. (2016) demonstrated that the relationship between the ELA (for which SLA is a surrogate) and the mass balance of Brewster Glacier is expected to be linear until ELA reaches the foot of the upper steep slopes at about 2000 m.The relationship then exhibits a marked inflection point indicating that changes in ELA become much more sensitive to changes in mass balance.Based on this result, the empirical relationship between SLA and extended mass balance was modelled by a piecewise linear regression to accommodate the change in sensitivity of the relationship at the foot of the upper steep slopes.In order to estimate uncertainties for the reconstructed mass balance while accounting for error in variables, unified linear regression was used (York et al., 2004), assuming an uncertainty of 30 m for the EOSS departure data.Departures for the full record were then substituted into the piecewise linear model to obtain estimates of annual mass balance over the EOSS observation period.These annual balance estimates (B EOSS a yr ) were then summed to obtain the cumulative mass balance of Brewster Glacier (1978-2013).The uncertainty envelope associated with the cumulative mass balance was obtained by a cumulative quadratic sum of the standard error associated with the estimated mass balance for each year as (3) 4 Results

Assessment of α MOD
Polar scatter plots in Fig. 7a, b illustrate the error α as a function of illumination and viewing geometry, using the isotropic and anisotropic correction models respectively.This does not reveal any substantial relationship to the solar zenith θ s .Nonetheless, it shows that the sensor view zenith θ v has a major influence on the accuracy of α MOD in both cases.Although the MODImLab albedo retrieval performs well at relatively low view zenith (θ v < 30 • ), its accuracy degrades markedly as θ v increases, in particular beyond 45 • .The use of the anisotropic reflectance model for snow and ice provides some, albeit limited, improvement (see Table 1), but does not mitigate the increasing underestimation of α MOD as θ v increases.These results are to be interpreted along with the fact that the AWS is located within a single 250 m pixel from the glacier margins (Fig. 5).It suggests that the influence of θ v on the degradation of performances is likely to be caused by the sensor panoramic distortion rather than a limitation of the albedo retrieval algorithm.As θ v increases, the pixel collocated with the AWS has a radiometry governed by a growing footprint that quickly includes darker rocks near the glacier margins (Figs. 1, 5).This causes spectral mixtures with a lower apparent reflectance resulting in a lower α MOD , explaining in turn the underestimation.This source of error can compromise the retrieval of albedo of pixels close to the glacier margins.In turn, the monitoring of glacier-wide albedo over a period of time during which glacier margins have changed may exhibit a trend that could jeopardize the mass balance model.Dumont et al. (2012) proposed adjusting the pixel radiometry using the outcome of the spectral unmixing (Sirguey et al., 2009) to ensure that the albedo retrieval only accounts for the relative share of radiation reflected by the snow and ice fraction.Despite its merit, this potential improvement was not tested here.Instead, the limitation of the albedo retrieval method was assumed to be mitigated by averaging values from pixels selected to be consistently and near-fully covered by snow and ice, as well as by considering a period of time during which the glacier margins have not substantially changed.
To appreciate better the performance of the albedo retrieval without this effect, Fig. 7c, d illustrates the agreement between α MOD and α AWS considering only θ v < 30 • , with the isotropic and anisotropic model respectively.Considering only clear-sky observations from both the MODIS cloud mask and AWS record (i.e.solid points in Fig. 7), Fig. 7c, d demonstrates that the albedo retrieval performs very well with no or negligible bias and a RMSD = 0.07 consistent with Dumont et al. (2012).The magnitude of the decrease of RMSD from 0.073 to 0.070 (Table 1) when using the anisotropic reflectance model is limited but compatible with results from Dumont et al. (2011Dumont et al. ( , 2012)).The better performance of the anisotropic reflectance model is strengthened Table 1.Agreement between α MOD and α AWS under various scenarios, namely with and without the anisotropic reflectance model for snow and ice, using MODImLab cloud mask vs. confident clear-sky days determined from incoming shortwave from the AWS, and considering view zenith up to 30 and 45 by the better fit to the 1 : 1 ratio line despite the apparent, albeit insignificant, decrease in R 2 from 89 to 86 % (p = 0.48).It is possible that the limited benefit of the anisotropic correction be due to the surface roughness of the surface being greater than that of the samples used to infer the anisotropic correction (Dumont et al., 2011) and possibly smoothing the anisotropic response of snow (Hudson et al., 2006).
Figure 7 shows that relying on the MODImLab cloud mask (i.e.empty and solid points) involves a number of larger errors that reduces the apparent accuracy of the method with RMSD increasing from 0.070 to 0.115 (Table 1).Such errors are, however, the result of less than 15 % MODIS observations with albedo estimates being grossly erroneous due to omission errors in the cloud mask.In view of this, we considered that the rate of good data points when relying solely on the MODImLab cloud mask was sufficient to appropriately capture the variability of glacier albedo despite the noise introduced by erroneous data.Such noise was mitigated by filtering the final sequence of glacier-wide albedo using a rolling average (see Fig. 8).Despite the reduced agreement between α MOD and α AWS when θ v up to 45 • are considered, the effect of growing pixel footprints with increasing view zenith is expected to be less pronounced on the glacier-wide albedo value α due to the widening of the glacier and the averaging effect.The compromise of retaining imagery with θ v < 45 • had the additional advantage of nearly doubling the number of points in the α time series compared to θ v < 30 • .

MODIS-derived albedo time series
Figure 8 illustrates the time series of retrieved glacier-wide albedo α(t), its yearly minimum during ablation season yr, namely α min yr , and the cumulative winter albedo, A w yr .From the time series of α(t), a mean annual cycle of glacierwide albedo for Brewster Glacier was obtained by averaging all values observed within each week of the year (Fig. 9).The standard deviation associated with the weekly distribution of albedo and the standard error of the corresponding mean (i.e. the standard deviation of the sample-mean estimate of this distribution) revealed the variability of the typical glacier-wide albedo cycle.Glacier-wide albedo of Brewster Glacier exhibits a pronounced seasonal cycle, with the highest albedos (greater than 0.7) observed between May and October.Glacier-wide albedo rises sharply from February to June until it reaches a maximum that coincides with the winter solstice (20/21 June).Air temperatures are lowest from June through to September; consequently fresh snow on the glacier surface remains in a less transformed state for longer.Albedo gradually decreases from September as a result of the metamorphism of snow and accumulation of light absorbing impurities.The decrease accelerates from Octo-ber to November as the glacier surface transitions from snow covered to predominately exposed ice (glacier-wide albedo ranging from 0.58 to 0.65).The glacier-wide surface albedo stabilizes from November through to December before another sharp decrease towards its minimum in February.
Initially it was expected that glacier-wide albedo would continue to decline through January and February, reaching a minimum at the end of the main melt season (late February/March).Instead, glacier-wide albedo on Brewster Glacier is shown to reach its minimum typically during the first 2 weeks of February (Fig. 9, Table 2).This suggests that transient snowfall and early accumulation in the upper region of the glacier is not uncommon from late February and may prevent the glacier-wide albedo from reaching a minimum later at the end of summer.Table 2 lists the timing and magnitude of the minima reached during each mass balance year.Ten out of 14 minima (70 %) were observed between mid-January and mid-February.Figure 8 shows that the two occurrences in early December 2002 and 2003 may be considered misleading because the albedo stabilized at a relatively high value throughout the summer, making the minima sensitive to the smoothing and inherent noise of the albedo time series.The values and timing of the secondary minima shown in Table 2 show that a slightly higher value was reached later in February in both years, and the timing was consistent with all other occurrences.Nonetheless, the primary values were retained for the rest of the analysis.
The standard deviation of weekly albedo values and their corresponding means (light and dark shading in Fig. 9 respectively) reveal that the average albedo of the glacier surface exhibits relatively little variation in winter.This indicates that, during the accumulation season when the glacier is fully snow covered, the average albedo of fresh snow and snow in early stages of metamorphism exhibit limited annual change in optical characteristics.Conversely, the glacierwide albedo is highly variable over the ablation season.This  is due to the higher variability in the timing and extent of the exposure of ice over the summer.

MODIS-derived mass balance models
We assessed the relationship between MODIS-derived albedo values (α min yr and A w ) and observed winter, summer, and annual balances (Fig. 10).Glacier-wide albedo minima on Brewster Glacier were found to correlate very well with the observed annual balances (R 2 = 0.93).This result is similar to that found by Dumont et al. (2012) on Saint-Sorlin glacier, and confirms the validity of using α min yr as a predictor for the mass balance of Brewster Glacier and possibly other glaciers in the Southern Alps of New Zealand.Annual albedo minima were also well correlated to summer bal-ance (R 2 = 0.83), although a better model (R 2 = 0.86) with half the standard error of prediction was found between the winter cumulative albedo and winter balance.
Figure 11a illustrates the agreement between MODISderived annual balances B MODIS a yr obtained with Model M1 and the observed annual balance from the glaciological method.The 93 % determination and a linear regression very close to (and not significantly different from) the 1 : 1 ratio line confirms that the record of mass balance for Brewster can be reliably extended over the MODIS observation period (i.e. 2000-2013).The standard error of prediction of B MODIS a yr of about ±0.5 m w.e.compares to that obtained from the glaciological method (±0.3 m w.e., see Cullen et al., 2016).Figure 11b shows that the estimation of BM12 s yr = BM1 a yr − BM2 w yr−1 also fits the observations very well, and better than if determined directly by α min yr (R 2 = 88 % vs. 82 %).This confirms that the variability of α min yr captures characteristics associated with both accumulation and ablation rather than ablation alone.This can be explained by the fact that α min yr is a proxy for the AAR, which is also controlled by winter balance.Alternatively, the cumulative winter albedo can be viewed as independent from ablation for it is measured in advance of ablation and involves processes that only favour accumulation.
The extended seasonal and annual balance record for Brewster Glacier based on MODIS-derived albedo data is presented in Fig. 12   = 92 %, p = 4.9 × 10 −8 ), thus confirming the findings of Cullen et al. (2016) who showed that 87 % of annual balance was controlled by summer balance.This finding compares favourably to observations in the French Alps where summer balance was also shown to be largely responsible for changes in the annual balance (Thibert et al., 2013;Six and Vincent, 2014).

Mass balance reconstruction with EOSS data
Brewster Glacier has an SLA record from the EOSS survey programme dating back to 1977.We compared the annual mass balance record extended over the period 2000-2013 to the EOSS record to assess whether variability in mass balance matched the variability of SLA departures from ELA 0 .The EOSS estimates, obtained from aerial photographs, are completely independent from the observed MODIS-derived mass balance measurements.
Figure 12 shows that the variability in reconstructed annual balance closely follows the variability in the end of summer snowline elevation despite a visible outlier in 2012.This outlier is likely an artefact of the approach used to assign SLA values from the oblique photographs.Over the 14year period, the snowline from the EOSS survey was only digitized on four occasions (i.e. 2000, 2007, 2009, 2010), and not always from rectified photos.For all other years, when the quality and/or timing of the photo was inadequate, the snowline altitude was estimated based on a subjective placement of the photos in the sequence and using values of years which bracket the year's snowline under consideration (Willsman et al., 2015).Using this method, the years 2008, 2011, and 2012 were assessed as being nearly identical by the EOSS survey.Nevertheless, Cullen et al. (2016) questioned the validity of the 2012 ELA estimate from the EOSS programme with regard to the glaciological record.The 2012 EOSS photograph exhibits an exposed glacier with a snowline visible at the bottom of the north-eastern slopes rather than close to the ridge as interpreted (i.e.2270 m a.s.l.).For the portion of the piecewise regression above 2000 m a.s.l., replacing the interpreted 2012 SLA value by that obtained from the regression against the ELA Alps (i.e.2127 m a.s.l., see Sect.2.4) significantly improved the correlation between SLA departures from ELA 0 and the extended mass balance   Glacier (2000Glacier ( -2013)).Note that the 2012 SLA appears twice; the empty point is the original departure estimate (335 m) while the solid point is that obtained from the relationship with ELA Alps (197 m).record from R 2 = 0.27 (p = 0.233) to R 2 = 0.55 (p = 0.057).The resulting piecewise linear model shown in Fig. 13 was used to reconstruct the mass balance record over the duration of the EOSS programme .The unified linear regression yielded a standard error of 0.48 and 0.39 m w.e. for the regression below and above 2000 m a.s.l.respectively, which is believed to be reasonably conservative for the reconstructed values.
The reconstructed 37-year-long mass balance record and corresponding cumulative mass balance are shown in Fig. 14.The 37-year reconstructed mass balance record documents, for the first time, a relatively long mass balance sequence of a New Zealand's glacier and sheds new light on its variability.Overall, it is marked by three dominantly positive sequences, namely 1977-1989, 1991/1998, and 2001/2007, interrupted by an outstandingly negative glacier balance recorded in 1990 and 2 consecutive years of significantly negative balances in 1999 and 2000.Over the first 30 years of the record, it is estimated that the glacier gained an average thickness of up to 14.5 ± 2.7 m w.e.This was offset by a marked shift towards negative balances, which have dominated since 2008.Over the 2008-2013 period only, the glacier lost 5.1 ± 1.2 m w.e., or 35 % of the gain accumulated over the previous 30 years, bringing the overall estimated gain from 1977 to 2013 to 9.4 ± 3.0 m w.e.values (1977-2013).Also shown is the cumulative mass balance of Brewster Glacier (1977Glacier ( -2013)), based on the annual records reconstructed from EOSS values.

Brewster mass balance record
The estimate of cumulative mass balance for Brewster Glacier should be treated with caution as there is currently no way to assess any systematic bias in the glaciological method that would propagate into the mass balance modelled using MODIS observations and, in turn, the reconstruction of mass balance via EOSS data.At this stage, the validation and possibly correction of the mass balance series using geodetic estimates of volume change is crucial and could follow the framework outlined by Zemp et al. (2013).
To make a first quantitative estimate of volume and mass change, a GPS survey conducted in January 1997 up to the foot of the steep north-eastern slopes (Willis et al., 2009) is compared to the 1986 topography derived from photogrammetry of imagery captured in February 1986 (see Columbus et al., 2011), which reveals substantial changes in surface elevation (Fig. 15).Over the 11-year period, the area above the 1730 m contour has exhibited a marked gain of up to +25 m, while the ablation area below the 1730 m contour appears to have down-wasted in excess of −25 m.Integrated across the surveyed area, the average change in volume is estimated to be +7.63 m, or +6.88 m w.e. using a density of ice at 900 kg m −3 .Over the same period, the cumulative mass gain from the EOSS model indicates +6.35 ± 1.64 m w.e.Despite the lack of information in the upper accumulation area, the 1997 survey supports the reconstructed mass gain of Brewster Glacier over this period, with such gain even possibly underestimated.More accurate geodetic estimates of volume/mass change should, however, be completed, for example based on archive stereo imagery from airborne (aerial survey on 13 March 2006) and very high-resolution spaceborne sensors.
The interpretation of oblique photographs from the EOSS programme provides some qualitative insight about the relevance of the cumulated balance estimates.From the first photo in 1978 to about a decade later in 1987, the terminus retreated rapidly and a proglacial lake formed, events which could question the accrued mass estimate.However, close observation of the middle and upper part of the glacier shows evidence of substantial mass gain, as ice margins appear to have risen significantly, even burying some rock outcrops in the steeper slopes.These observations are, therefore, congruent with the surface elevation changes presented above.
The fast terminus retreat observed during the first 10 years of the EOSS survey photos supports the hypothesis that a dominant sequence of negative mass balance years occurred pre-1978, with retreat being the delayed response to these previous negative mass balances.A shift from negative to positive mass balance around 1978 is coincident with a major shift in the Pacific (Inter) Decadal Oscillation (PDO/IPO) from a cool to a warm phase in 1977/1978 (Salinger et al., 2001), incidentally also the start of the EOSS programme.Figure 16 shows that the series of PDO/IPO since 1978 corresponds well with the reconstructed mass balance of Brewster Glacier.The dominantly neutral to positive mass balance of Brewster Glacier between 1978 and 1998 is consistent with the warm phase of PDO/IPO, as well as with observations of advance/positive mass balance sustained by other glaciers in the Southern Alps (Chinn et al., 2005).The reversal to a dominantly cool phase after 1998 corresponds to the observed shift towards more frequent moderate to strong negative balances, interrupted by gains over the 2003-2007 period, during which the PDO/IPO also returned to a neutral or warmer phase.The correspondence between mass balance of Brewster Glacier and the phase of the PDO/IPO suggests that observed mass gain of Brewster Glacier since 1977 is only offsetting a longer sequence of dominantly negative balances  Despite the compelling correspondence of decadal trends in mass balance and the PDO/IPO, the mechanisms controlling mass balance over decadal timescales are still unclear.The connection between large-scale modes of the oceanatmosphere system and mass balance in the Southern Alps has been explored by numerous authors (e.g.Fitzharris et al., 1997Fitzharris et al., , 2007;;Tyson et al., 1997;Clare et al., 2002;Purdie et al., 2011), revealing the strength and position of the westerly circulation to be a key control.However, correlations between individual climate variables (e.g.air temperature and precipitation) and mass balance anomalies do not provide a full explanation of mass balance variability (e.g.Purdie et al., 2011).Further work is needed to explicitly resolve the mechanisms responsible for driving mass balance anomalies through modelling of the surface energy and mass balance of Brewster Glacier over the EOSS period.Given the dominant role of precipitation phase -albedo feedbacks in controlling mass balance sensitivity to air temperature (Conway and Cullen, 2016), careful consideration needs to be given to changes in the amount and timing of snowfall and the magnitude of the associated feedbacks.In addition, trends in cloud cover should be examined, as cloud cover has also been shown to significantly alter the relationship between air temperature and melt (Conway and Cullen, 2016).

P. Sirguey et al.: Reconstructing the mass balance of Brewster glacier
In the later part of the EOSS photo series from 2008 to 2015, there is evidence of substantial lowering in the upper part of the glacier, thus giving visual evidence of the recent sequence of negative balances.As terminus retreat has continued unabated over the time period of sustained significant mass gain, it unlikely that an advance can be anticipated in response to the mass gain.A shift in climate and enhanced ablation in the lower part may have hampered the ability of Brewster Glacier to recover area via ice transport.If future negative mass balances recur frequently, through either a return to a cool phase of the PDO/IPO or regional anthropogenic climate forcing, we may further anticipate a possibly accelerated retreat of the terminus in the near future, and a significant readjustment of the geometry of Brewster Glacier to the new climate.The application of a glacier flow model would allow further exploration of the climate/ice dynamics and/or glacier geometry relationships of Brewster Glacier.

On the suitability of α min yr in the New Zealand context
The successful application of MODIS-derived minimum glacier wide albedo (α min yr ) as a proxy for mass balance on Brewster Glacier indicates the potential for this method to be extended at a broader scale within the Southern Alps.Long-term records (2000-present) of α min yr variability could be generated for any glacier in the Southern Alps, providing that they are of a sufficient size (> 2 km 2 ) not to compromise the albedo retrieval from the relatively coarse MODIS pixels.Even if the glaciers do not have observed mass balance data from which to derive quantitative mass balance values, there is merit in resolving long-term trends.This would provide a method of tracking mass balance trends that is independent from the EOSS programme, with the benefit of less subjectivity and fewer logistical requirements.Trends in α min yr could be used to establish equilibrium status of a large number of glaciers, and assess whether negative mass balances observed at the 50 index glaciers extend to any of the other 3000 glaciers in the Southern Alps.Similarly, the timing and magnitude of α min yr on glaciers distributed throughout the Southern Alps could reveal new insights about spatial contrasts in the behaviour of New Zealand glaciers, as well as about their relationship to local and global climate.
At present, only 50 glaciers (out of the more than 3000 inventoried in New Zealand) are the subject of EOSS/SLA monitoring with a flight date being timed for the first suitable flying weather after 1 March (Chinn et al., 2005).However, the seasonal variability of glacier-wide albedo on Brewster Glacier revealed that the expected date for the minimum to be reached is within the first 2 weeks of February.Although albedo may rise due to new snow at high elevation, not necessarily concealing the snowline, this suggests that new snow as early as the second half of February is not exceptional on Brewster Glacier.Exploring the spatial distribution of this timing could inform the EOSS programme about the appropriate scheduling of the survey to maximize chances of resolving the end-of-summer snowline.
The changing albedo of glaciers and ice sheets is currently a topic of huge interest, with a particular focus on the Greenland Ice Sheet (Mernild et al., 2014;Dumont et al., 2014) and western European glaciers (Gabbi et al., 2015;Naegeli et al., 2015).Although the present study did not explore trends in the albedo sequence in detail (e.g.albedo of winter snow only), by being situated in the path of the prevailing westerly airflow, the Southern Alps of New Zealand are exposed to deposition of ashes and dust transported from Australia across the Tasman Sea.As the MODIS albedo retrieval allows additional optical properties of snow and ice to be estimated, such as specific surface area (SSA) (Mary et al., 2013) and effective impurity content (e.g.Black Carbon and mineral dust), there is potential to explore further the albedo sequence on selected New Zealand glaciers and snowfields.Trends in albedo and controlling factors such as SSA or impurities could shed light on processes affecting snow in New Zealand and their possible feedback on glacier variability.

On the estimation of winter accumulation from albedo
This study successfully tested an approach to estimate winter accumulation on a glacier by tracking the accumulated winter albedo.Although the value of the optimum albedo threshold found here (0.72) needs to be confirmed on other glaciers, it shows that the glacier-wide albedo signal retrieved from MODIS imagery holds useful information about the variability in accumulation.The metric may also be useful in the study of snowpacks outside glacierized areas, perhaps providing real-time information about the status of seasonal snow.

Conclusion
In this study, we assessed the validity of the MODImLab albedo retrieval algorithm using the 2-year-long record of albedo on Brewster Glacier, New Zealand.As part of the accuracy assessment, we tested the sensitivity of MODISderived albedo measurements to illumination/viewing geometry, anisotropic reflectance of snow and ice, and omission errors in the cloud masking algorithm.The results show that the main source of error originates from the growing footprint of MODIS pixels with increasing viewing angles that yield an underestimation of retrieved albedo, likely due to spectral mixture with nearby rocks.Overall, the accuracy of the albedo retrieval was found to decrease from RMSD = 0.07 in clear-sky conditions to RMSD = 0.115 when relying on the MODImLab cloud mask.Despite gross errors due to cloud omission, the reliance on the MODImLab cloud mask did not appear to compromise the albedo sequence, confirming it as a useful method to retrieve values of glacier-wide albedo at a reasonable frequency through a season.MODIS-derived albedo minima correlated well (R 2 = 0.93) to the observed annual balance of Brewster Glacier, confirming that minimum glacier-wide albedo can be an appropriate proxy for mass balance in this maritime environment.Further to using MODIS-derived albedo minima to extend the mass balance record of Brewster Glacier over the period of MODIS observations (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013), we proposed an alternative mass balance model that enabled the reconstruction of separate winter and summer balance records.In particular, cumulative winter albedo, defined as the integral of the glacier-wide albedo during winter when such albedo exceeds a given threshold, proved to correlate strongly with winter balance (R 2 = 0.88).Although the glacier-wide minimum albedo also correlated well to the summer balance (R 2 = 0.82), the latter was best estimated from annual and winter balance (R 2 = 0.88).This is an important development of the albedo method to estimate seasonal glacier mass balance from multi-spectral imagery as the winter and summer reconstruction approach has the potential to provide detailed information on accumulation and ablation processes.Although the albedo method to measure mass balance still requires glaciological measurements for calibration, this confirms the validity of monitoring surface albedo of mountain glaciers to retrieve mass balance signals.This is promising as it provides the platform from which to assess the behaviour of other glaciers with relatively little logistical cost.The extended record also confirms that annual balance of Brewster Glacier is largely controlled by the summer balance as reported by Cullen et al. (2016)  (R 2  = 92 %, p = 4.9 × 10 −8 ).The strong correlation between the extended annual balance record and Brewster Glacier's EOSS/SLA record allowed the annual balance to be reconstructed over the period 1977-2013, thus providing the longest mass balance record for a glacier in New Zealand.Between 1977 and1997, Brewster Glacier exhibited dominantly positive mass balances only interrupted in 1990 by an outstandingly negative balance.This period of positive mass balance is consistent with the behaviour of other glaciers in the Southern Alps and stresses the high sensitivity of New Zealand glaciers to shifts in regional ocean-atmospheric state (Fitzharris et al., 2007).Over the 37-year period, our results show that Brewster Glacier gained 14.5 ± 2.7 m w.e. of mass between 1978 and 2007.Although estimates of the accrued mass could not be validated due to the lack of comprehensive geodetic estimates, a GPS survey in 1997 over most of the glacier compared to the 1986 topography confirmed the magnitude of the gain over a 11-year period.This gain was offset by a marked shift toward negative balances after 2008, yielding a loss of 5.1 ± 1.2 m w.e., or 35 % of its gain accumulated over the previous 30 years.The future evolution of Brewster Glacier will ultimately depend on how anthropogenic and internal forcing of local and regional climates are transferred to sur-face energy and mass balance anomalies and how these will influence the dynamic response.Further work to resolve the mechanism for driving the mass balance anomalies is still required, ideally by modelling the surface energy and mass balance of Brewster Glacier over the EOSS period.

Figure 6 .
Figure6.Measured SW↓ compared to that predicted by integrating the SPCTRAL2 clear-sky solar spectral model ofBird and Riordan (1986) over the range of wavelengths captured by the pyranometer i.e. (300-2800 nm) over 4 days in summer (a) and winter (b).Solid triangles correspond to MODIS observations retained in the strict subset of images with guaranteed clear-sky conditions.Empty triangles would be discarded from this subset even if the MODImLab cloud mask had failed to detect clouds.

Figure 7 .
Figure 7. Agreement between MODIS-derived albedo (α MOD ) and observed albedo (α AWS ) using the isotropic (left) and anisotropic (right) reflectance model.Polar scatter plot of all clear-sky MODIS observations over the AWS measurement period showing the quality of the albedo retrieval as a function of solar and view zenith (top).The black symbols correspond to points with highly confident clear sky determined from the AWS record.Agreement between α MOD and α AWS for points with view zenith θ v < 30 • (bottom).The dotted line is the ideal 1 : 1 ratio.The dashed lines show the confidence interval for the linear regression (1σ ).

Figure 9 .
Figure 9.The average seasonal albedo cycle for the study period (2000-2013) presented as a mean weekly glacier-wide albedo.The dates of the annual albedo minima are indicated by the vertical lines.The darker and lighter envelopes around the mean indicate the standard error and standard deviation respectively of the set of observations used to compute the mean.

Figure 10 .
Figure 10.Unified linear regression between MODIS-derived minimum albedo and (a) annual balance, (b) summer balance, along with (c) the relationship between cumulative winter albedo and winter balance.Note the quality of the linear regression models conveyed by the coefficient of determination R 2 and the standard error of prediction.The dashed line is the standard error prediction interval.

Figure 11 .
Figure 11.Relationship between MODIS-derived balance for (a) annual balance with Model M1 and (b) summer balance with Model M12.The dotted line is the 1 : 1 ratio.The dashed line is the standard error prediction interval.

Figure 12 .
Figure 12.MODIS-derived seasonal and annual mass balance series of Brewster Glacier (solid line) compared to estimations from glaciological observations (dotted line).Standard errors of prediction are shown in shaded envelopes.Note the record of end-ofsummer snowline elevation (SLA) from the EOSS programme and its close match to the annual balance series.

Figure 13 .
Figure13.MODIS-derived annual balance vs. EOSS departures from ELA 0 for BrewsterGlacier (2000Glacier ( -2013)).Note that the 2012 SLA appears twice; the empty point is the original departure estimate (335 m) while the solid point is that obtained from the relationship with ELA Alps (197 m).

Figure 14 .
Figure 14.Comparison between the observed annual balance measurements (2005-2013) and the two reconstructed annual balance records.The first reconstruction (Model M1) is based on minimum albedo values (2000-2013), while the second reconstruction is based on EOSSvalues (1977-2013).Also shown is the cumulative mass balance of BrewsterGlacier (1977Glacier ( -2013)), based on the annual records reconstructed from EOSS values.

Figure 15 .
Figure 15.(a) Map of change in elevation between interpolated GPS transects captured in January 1997 (Willis et al., 2009) and 15 m DEM from large-scale national topographic mapping based on imagery dated February 1986 (Columbus et al., 2011).(b) Change in profile between 1986 and 1997 ( H are exaggerated 5 times).
s yr , where BM2 s yr = s 2 α min yr + o 2 and BM2 w yr = s 3 A w yr + o 3 , Sirguey et al.: Reconstructing the mass balance of Brewster glacier of seasonal mass balance consistent with BM1 a yr , B w was obtained from Model M2 while B s was deducted as BM12 s yr = BM1 a yr − BM2 w yr−1 .The standard errors associated with BM12 s yr values were obtained via Gaussian error propagation.Since MODIS winter observations only started in 2000, Model M2 could only support estimates of BM2 www.the-cryosphere.net/10/2465/2016/The Cryosphere, 10, 2465-2484, 2016 P.
Figure 8. Time-series of MODIS-derived glacier-wide albedo α (t).The shaded sections indicate the areas under which the signal was integrated to measure the cumulative winter albedo value A yr .The contrasting annual minimum albedo values α min yr are also evident.

Table 2 .
Dates and minima of the annual glacier-wide albedo.The dates and values of the secondary minimum for mass balance years 2002/2003 and 2003/2004 are shown in italics.