Applying a distributed mass-balance model to identify uncertainties in glaciological mass balance on Brewster Glacier, New Zealand

Abstract

A distributed mass-balance model is used over a 10-year period for the re-analysis of a glaciological mass-balance time series obtained from Brewster Glacier, New Zealand. Mass-balance modelling reveals glaciological mass balance has been overestimated, with an average mass loss of −516 mm w.e. a⁻¹ not captured by observations at the end of the ablation season, which represents 35% of the annual mass balance. While the average length of the accumulation season (199 days) remains longer than the ablation season (166 days), melting is not uncommon in the core part of the accumulation season, with 2–32% of total snowfall being melted. Refreezing of meltwater is also important, with 10% of surface and subsurface melt being refrozen in the present climate. Net radiation, driven primarily by net shortwave radiation, is the main contributor to melt energy, with melt variability mainly influenced by the turbulent heat fluxes, net longwave radiation and the heat flux from precipitation in the ablation season. Snowfalls in summer are an important moderator of melt, highlighting the critical role of the ice-albedo feedback and phase of precipitation on seasonal mass balance. A complete homogenisation of the long-term glaciological mass balance for Brewster Glacier is still required.

1. Introduction

Glaciers in the Southern Alps of New Zealand have undergone significant changes over the past few decades, which have been observed in both remote-sensing and field measurements (e.g. Sirguey and others, 2016; Cullen and others, 2017; Hugonnet and others, 2021; Lorrey and others, 2022). Surface thinning rates have increased almost seven times over the period 2015–2019 compared to 2000–2004, with a record value of 1.52 ± 0.5 m a⁻¹ (Hugonnet and others, 2021). Accelerated mass loss of the Southern Alps glaciers has been linked to atmospheric temperature changes, which will likely become more evident as global warming continues (Vargo and others, 2020; Lorrey and others, 2022). Although the acceleration of glacier mass loss shows connections to a warming atmosphere, decadal rates of glacier mass change remain strongly modulated by regional climatic conditions (Hugonnet and others, 2021). For example, during the period of 1983–2008, 58 glaciers in the Southern Alps were believed to advance while glaciers globally retreated (Mackintosh and others, 2017). Mackintosh and others (2017) argued that this phase of glacier advance, however, was mainly due to reduced regional air temperatures associated with southerly winds and low sea surface temperatures in the Tasman Sea. The contrasting patterns of global glacier behaviour continues to be a major topic of research, with emphasis on regional climatic conditions and the relative importance of temperature and precipitation controls (Mackintosh and others, 2017; Wang and others, 2019; Hugonnet and others, 2021; Tielidze and others, 2022).

Glacier mass balance is primarily controlled by the energy and mass exchanges between the atmosphere and the glacier surface. Due to their high sensitivity, glaciers are recognised as direct indicators of climate change, with any such change being observed in the annual glacier mass balance (Chinn, 1996). The implications of glacier fluctuations to the hydrological characteristics of glacierised drainage basins, coupled with the contributions to global sea level rise, led to a global initiative of the monitoring of glacier mass balance. Standardised observations of global glacier changes are collected and compiled by the World Glacier Monitoring Service (Müller and others, 1977; Haeberli, 1998; Zemp and others, 2009). Global estimates of glacier mass balance are primarily derived from the glaciological method, where point measurements of accumulation and ablation across the glacier are obtained (Østrem and Brugman, 1991). While there are several methods to assess glacier mass balance, the glaciological method is one of the most widely used and is referred to as the direct method (Pelto, 1988; Valla, 1989; Kuhn and others, 1999; Sicart and others, 2007; Cullen and others, 2017). This method provides an important baseline for the monitoring of glacier mass balance but is often subject to systematic and random uncertainties resulting in small annual but large cumulative errors (Andreassen and others, 2016; Wagnon and others, 2021). To address the uncertainties, the glaciological dataset is sometimes calibrated and combined with geodetic volume changes (Cogley, 2009; Huss and others, 2009; Thibert and Vincent, 2009; Zemp and others, 2009; Marinsek and Ermolin, 2015; Thomson and others, 2017) and numerical model outputs (Jóhannesson and others, 1995; Schuler and others, 2005; Mölg and others, 2008). While the glaciological method provides seasonal to annual point-based measurement of glacier surface mass balance, the geodetic approach yields greater spatial coverage of surface elevational changes (usually surveyed over a few years to a decade), and is sometimes used to estimate the surface, basal and internal mass balance of glaciers (Andreassen and others, 2016). The geodetic method is particularly useful in detecting biases in glaciological measurement campaigns and calibrating cumulative glaciological mass balances (Cox and March, 2004). Furthermore, while the glaciological method provides the seasonal and annual surface mass balance, the time period associated with these mass-balance quantities are important – for example, any amount of melt after the end of summer survey is not included in the glaciological mass balance (Huss and others, 2009). Since the geodetic method does not provide estimates of interannual variability, numerical models (including simple temperature index models and physically based surface energy and mass-balance models) are valuable tools that are often used to simulate the seasonal and annual glacier mass balances and to test the sensitivity of glaciers to climate change (Oerlemans and Hoogendoorn, 1989; Oerlemans, 1997; Anderson and others, 2010; Arndt and others, 2021; Kinnard and others, 2022). However, numerical models rely on high-quality climate data, which can be challenging to obtain in some glaciated regions. The uncertainties and limitations inherent to each method suggest that comprehensive monitoring strategies include a combination of methods, such as glaciological observations, periodic geodetic surveys and modelled mass-balance estimates (Machguth and others, 2006; Mölg and others, 2008; Huss and others, 2009; Thibert and Vincent, 2009; Elagina and others, 2021; Kinnard and others, 2022).

In the Southern Alps of New Zealand, monitoring glacier behaviour has been done using several approaches, including the end of summer snowline (EOSS), glaciological method, numerical modelling and remote-sensing techniques. Brewster Glacier, located just west of the main divide in the Southern Alps (Fig. 1), is one of the index glaciers in the EOSS survey. The EOSS survey captures low altitude oblique aerial photographs from aircraft for more than 50 index glaciers, providing insights about interannual variations in snowline over glaciers in the Southern Alps (Willsman and others, 2015; Lorrey and others, 2022). Brewster Glacier is a benchmark glacier for research (e.g. Smart and others, 2000; Chinn and others, 2005; George, 2005; Anderson and others, 2010; Gillett and Cullen, 2011; Conway and Cullen, 2013, 2016; Conway and others, 2015; Cullen and Conway, 2015; Vargo and others, 2020; Abrahim and others, 2022) and holds the longest in situ record of glaciological mass balance in the Southern Alps (Sirguey and others, 2016; Cullen and others, 2017). Cullen and others (2017) re-analysed the glaciological observations and developed a geostatistical method to calculate mass balance, which allowed the spatial and temporal variability in glacier mass balance of Brewster Glacier to be resolved. During 2005–2015, Brewster Glacier recorded an annual glacier-wide mass balance of −102 mm w.e., indicating that Brewster Glacier has been losing mass during this period (Cullen and others, 2017). The glaciological mass balance has provided a baseline to apply other methods to reconstruct mass balance, including the surface albedo method using moderate resolution imaging spectroradiometer (MODIS). Sirguey and others (2016) showed that minimum glacier-wide albedo obtained from MODIS imagery is a reliable predictor of summer mass balance and regular monitoring of glacier-wide albedo also makes it possible to determine winter accumulation, enabling the seasonal mass-balance record of Brewster Glacier to be extended back to the start of MODIS observations in 2000. To reconstruct mass balance even further back in time, Sirguey and others (2016) established a statistical relationship between the end of summer snowlines from the EOSS programme and the MODIS derived record of mass balance. The 37-year reconstructed mass balance for the period 1977–2013 shows a cumulative mass gain through to 2007, followed by a marked shift to mass loss after 2008. This latter period resulted in a 35% loss of the accumulated mass gained over the past 30 years (Sirguey and others, 2016).

Fig. 1. (a) An oblique view of the location of Brewster Glacier in the Southern Alps of New Zealand; and (b) map of Brewster Glacier showing the location of the AWS and surrounding topography. Contour lines are at 100 m intervals. The glacier margin is based on the 2011 survey. The probe locations are shown in the black dots while the stake locations are shown in the black circle. The red dots represent points that are no longer measured (4 November 2021). The insert map depicts the location of Brewster Glacier within New Zealand.

To revisit this assessment, Sirguey and others (2019) reanalysed historical aerial surveys with modern photogrammetric methods to obtain a geodetic assessment of mass loss over the periods 1986–2006 and 2006–2018. The revised record showed slow cumulative mass loss over the 1986–2006 period, followed by strong significant mass loss between 2006 and 2018. These observations immediately raised doubt as to whether the glaciological mass balance is capturing the full extent of the mass loss on Brewster Glacier, given the moderate cumulative mass loss shown by the direct method over the same period. Sirguey and others (2019) recalibrated the mass-balance record for Brewster Glacier and clearly demonstrated that the uncalibrated glaciological record underestimates mass loss, highlighting the need to reconsider uncertainties in the glaciological method. While previous studies have shown mass loss occurs after end of summer mass-balance observations are made (Cullen and Conway, 2015; Conway and Cullen, 2016), the full extent of glacier-wide mass loss for the entire ablation period has not yet been resolved.

To address uncertainties in the glaciological mass-balance record, this study uses a mass-balance model to assess how near-surface climate influences the spatial and temporal variability in energy and mass exchanges over Brewster Glacier. While Anderson and others (2010) provided useful insights on the glacier-wide mass balance of Brewster Glacier over a 4-year period, this is the first to use a fully distributed, physically based mass-balance model over an extended period on Brewster Glacier. This study aims to (i) reconstruct the mass balance of Brewster Glacier over a 10-year period using a calibrated physically based model, (ii) quantify the magnitude of mass loss that occurs outside the glaciological observation period, and (iii) detail the energy fluxes associated with periods of extreme mass loss. Meteorological data from an automatic weather station (AWS) and glaciological mass-balance measurements are used in the mass-balance modelling. The modelling framework takes advantage of the known meteorology of the study site to provide important details about the role each energy flux plays in influencing the variability in seasonal mass balance.

2. Study site and data

2.1. Site description

Brewster Glacier is a small mountain glacier in the Southern Alps, which in March 2011, covered an area of ~2.03 km² and descended from 2389 to 1706 m a.s.l. (Fig. 1). Brewster Glacier experiences a temperate, maritime climate characterised by high levels of precipitation (~6000 mm w.e.). Annual glacier-wide air temperature (T _a) is estimated to be equal to 0.6°C between 2010 and 2020. At least 77% of the total glacier area lies below 2000 m a.s.l. and has a slope of 10°, while the upper region is steeper with an average slope of 31°. Additional details of the study site are provided in Abrahim and others (2022). A 50 m digital elevation model (DEM) is used to represent the surface topography of Brewster Glacier and surrounding topography for the mass-balance modelling component of this research, which has been resampled from a 15 m resolution DEM from NZSoSDEM (Columbus and others, 2011). The NZSoSDEM is interpolated from 20 m contours from the Land Information New Zealand (LINZ) 1:50 000 topographic map derived from a photogrammetric aerial survey in February 1986 (Columbus and others, 2011). The glacier boundary is determined using a near snow-free QuickBird image from the KiwImage project obtained on 8 February 2011. Oblique aerial photographs taken on 13 and 30 March 2011 are also used to define the glacier boundary.

2.2. Meteorological data

A meteorological dataset obtained over the period June 2010 to November 2020 from an AWS located next to the terminus of Brewster Glacier at 1650 m a.s.l. is used in this study (previously referred to as AWS_lake in Cullen and Conway (2015) and Abrahim and others (2022); Fig. 1). A full description of the instruments along with accuracy, data ranges and associated uncertainties is provided in Abrahim and others (2022). It should be noted that due to sensor issues, wind speed data are unavailable for the period of 1 January 2018 to 1 November 2019. During this period, average wind speed at Brewster Glacier is used. Average annual air temperature recorded at AWS_lake is 2.1°C and the seasonal cycles of surface temperature (T _s) and albedo (a) indicate that the site is snow-covered between June and November. The atmosphere is humid and frequent movements of high- and low-pressure systems occur across the region. Average wind speed (U) is 3.3 m s⁻¹ and cloud cover is frequent with a distinct diurnal and seasonal cycle (Abrahim and others, 2022).

2.3. Glaciological data

Glaciological observations of the seasonal surface mass balance of Brewster Glacier commenced in February 2004. The first year of measurements are described by George (2005). The glaciological method is used to measure the summer and winter mass balances. At the end of winter (usually early to mid-November), a field campaign is conducted where geolocated point-based measurements of snow depth (using snow probes) are obtained (Fig. 1, probe locations shown in black dots), while snow density is measured from at least one snow pit that is dug on the centre line of the glacier. It should be noted that the probing locations above 2000 m a.s.l. are not regularly sampled due to the logistical challenges of accessing the steeper parts of the glacier. Some probing and stake locations are no longer measured as they now lie outside of the glacier margin (Fig. 1, red dots). The end of summer observations are usually obtained in March and include measurements of ablation from stakes drilled into the glacier surface during the end of winter field campaign (Fig. 1, stake locations shown in black circles). The glaciological mass balance is calculated using a floating-date system and therefore, the time period of each mass-balance year varies (Cullen and others, 2017).

To reduce spatial and temporal artefacts introduced by changes to the measurement regime as the mass-balance programme evolved, a geostatistical method is used to calculate mass balance (Cullen and others, 2017). The similarity in the spatial variability of the winter and summer balances over Brewster Glacier between an individual mass-balance year and other years allowed indices to be determined for both. Importantly, the indices are better predictors than using elevation in revealing the spatial structure of winter and summer balance. The geostatistical method has been used to update the glaciological mass balance to the present.

3. Methodology: mass-balance modelling

3.1. Mass-balance model

The mass-balance model developed by Mölg and others (2008) and already used on Brewster Glacier at a single point in the ablation zone (Conway and Cullen, 2016) is optimised to run in a fully distributed manner over the period 2010–2020. The mass balance (MB) is calculated as expressed in Eqn (1),

(1)$$MB = M_{{\rm acc}} + M_{{\rm ref}}-M_{{\rm melt}} + M_{{\rm dep}}-M_{{\rm sub}}-M_{{\rm subsurfmelt}} + {\rm d}W_{{\rm liq}}, \;$$

where M _acc is the snow accumulation, M _ref is the refrozen meltwater in the snowpack, M _melt is the surface melt, M _dep is the surface deposition, M _sub is the surface sublimation, M _subsurfmelt is the subsurface melt and dW _liq is the change in liquid water storage. Melt energy (Q _M), when T _S = 273.15 K at a 30 min time step, is calculated as,

(2)$$Q_{\rm M} = SW_\downarrow ( {1-a} ) + LW_\downarrow -\sigma \varepsilon T_{\rm S}^4 + Q_{\rm S} + Q_{\rm L} + Q_{{\rm PRC}} + Q_{\rm C} + Q_{{\rm PS}}, \;$$

where SW _↓ is the incoming shortwave radiation, a is the albedo, LW _↓ is the incoming longwave radiation, σ is the Stefan Boltzmann constant (5.67 × 10⁻⁸ W m⁻²), ɛ is the emissivity of snow/ice (assumed to be unity), Q _S is the sensible heat flux, Q _L is the latent heat flux, Q _PRC is the heat flux from rain, Q _C is the conductive heat flux and Q _PS is the penetrating shortwave radiation to the subsurface. In the model, energy fluxes directed towards the glacier surface are positive.

Input variables for the model include precipitation (P), T _a, relative humidity (RH), air pressure (p), U and cloud cover (N_ɛ), which are sourced from meteorological observations (Section 2.2.). These meteorological observations are collected from an AWS located on bedrock. The location of the AWS (not on the glacier surface) has implications on air temperature, particularly summer air temperatures (air temperature sensor height is 3.0 m), which have been corrected for solar heating following the approach used by Cullen and Conway (2015) and Abrahim and others (2022). Air temperature is modified when distributed across the glacier surface, while no vertical gradients are applied for wind speed and precipitation. The DEM discussed in Section 2.1 is used to spatially distribute air temperature using a lapse rate approach and to create a shading file for the estimation of incoming shortwave and longwave radiation across all elevations. Uncertainties relating to the parameterisation of surface roughness lengths (used in the exchange coefficient parameterisation for the calculation of the turbulent heat fluxes) and cloud effects on radiation and melt energy, as well as the calculation of surface temperature and albedo have been addressed in a series of point-based studies (e.g. Conway and Cullen, 2013, 2016; Conway and others, 2015; Cullen and Conway, 2015), which have helped to limit the uncertainty in the modelling framework used in this study.

The albedo parameterisation is based on three albedo factors including the albedo of ice (a _ice), firn (a _firn) and fresh snow (a _frsnow), where an e-folding constant (t*) is used to describe the change from fresh snow to firn (Oerlemans and Knap, 1998). The scheme has been slightly altered to ensure albedo returns to the albedo of the underlying snow or ice surface after fresh snowfall. The daily threshold for total snowfall depth has been modified to a value of 5 cm, which represents the threshold above which new snowfall affects albedo. Small quantities of snowfall are often redistributed into crevasses and hollows and will have a minor effect on the albedo (Conway and Cullen, 2016). Incoming shortwave and longwave radiation are calculated using cloud cover (N_ɛ), air temperature, air pressure and relative humidity (Conway and others, 2015). A full description of this method is presented in Abrahim and others (2022).

The turbulent heat fluxes (Q _S and Q _L) are calculated using the bulk aerodynamic method described in Conway and Cullen (2013). Turbulent sensible and latent heat fluxes are determined using:

(3)$$Q_{\rm S} = c_{\rm p}\rho _0\displaystyle{\,p \over {\,pa}}CU_Z\,( {T_Z-T_S} ) $$

(4)$$Q_{\rm L} = 0.622\;L_{\rm v}\rho _0\displaystyle{1 \over {\,pa}}CU_Z\,( {e_Z-e_0} ) , \;$$

where c _p is the specific heat capacity of air at constant pressure (1005 J kg^–1 K^–1), ρ ₀ is the density of air at standard sea level (1.29 kg m^–3 at 0°C), p is the air pressure, pa is the air pressure at standard sea level (1013 hPa), C is a dimensionless exchange coefficient, U_Z, T_Z and e_Z are the U (m s⁻¹), T _a (K) and vapour pressure (hPa) at height Z (m), T _S is the surface temperature (K), e ₀ is the vapour pressure at the surface (hPa) and L _v is the latent heat of vaporisation (2.514 MJ kg^–1), or the latent heat of sublimation (2.848 MJ kg^–1) if T _S is below the melting point.

Precipitation (P _scaled) is determined from the nearby Makarora AWS (30 km to the southwest of Brewster Glacier at an elevation of 320 m a.s.l.) using the scaling factors described in Cullen and Conway (2015). A tipping-bucket rain gauge is used to measure summer precipitation at AWS_lake. Due to snow-cover of the rain gauge, precipitation during winter and spring is not measured. A gauge undercatch factor of 1.25 is applied to the AWS_lake precipitation data. To develop a precipitation record for Brewster Glacier, Cullen and Conway (2015) fitted a linear least square relationship to precipitation at AWS_lake and Makarora AWS during the summer 2010/11. The precipitation record used in this study is based on this scaling factor. The partitioning of snow versus rain is determined using a snow/rain threshold range (e.g. Mölg and Scherer, 2012; Mölg and others, 2012) (−1 (below which all precipitation is solid), 1 (above which all precipitation is liquid)). The transition from solid to liquid is linearly interpolated and accounts for mixed phase precipitation. The fresh snow density is assumed to be 300 kg m⁻³ (Gillett and Cullen, 2011; Cullen and Conway, 2015). Heat from rain is also calculated using P _scaled (e.g. Bogan and others, 2003; Mölg and others, 2012). The iterative energy-balance closure method employed by Mölg and others (2008) was used to calculate T _S and is governed by the energy available at the glacier surface under the requirement of equilibrium among the surface energy fluxes. When surface temperature surpasses melting point, it is reset to 273.15 K and the resultant energy is the energy available for melt. Energy storage in the surface is included (Garratt, 1992) and T _S is set to 273.15 K to initialise the iterative scheme.

The subsurface and energy balance are connected through T _S, which determines the englacial temperature (T _en) on a numerical grid. The subsurface profile numerically resolves the thermodynamic energy equation through a 14-layer vertical subsurface module. In a layer with snow and ice, a weighted average using the values for the thermal diffusivity, thermal conductivity and density of snow and ice is used to resolve the thermodynamic energy equation. Penetrating shortwave radiation (Q _PS) and englacial refreezing processes are factors that are considered during the T _en calculations. The penetration of shortwave radiation into the subsurface is governed by constant fractions for snow and ice surfaces and is attenuated exponentially with depth (Table 1). For an ice surface, the fraction derived by Mölg and others (2008) is used, while the value for a snow surface is based on Bintanja and Van Den Broeke (1995). The extinction coefficients for snow and ice are based on Bintanja and Van Den Broeke (1995). While the subsurface routine has been used in previous studies (e.g. Mölg and others, 2012, 2014), these parameterisations have not been validated for Brewster Glacier and are subject to uncertainty. The conductive heat flux, which is the conductive heat flux between the surface and the upper layer of the 14-layer structure is determined using T _en. The subsurface temperature initialisation is isothermal with the initial profile set to 273.15 K. The bottom temperature of the subsurface structure is set as 273.15 K (e.g. Conway and Cullen, 2016); however, this assumption may contain some uncertainty particularly for the higher accumulation areas, where a 273.15 K bottom temperature can lead to an overestimation of melt, linked to lower air temperatures and snow cover (Pellicciotti and others, 2009). The density of the snowpack (settled snow = 500 kg m⁻³; ice = 900 kg m⁻³) is based on snow-pit measurements at Brewster Glacier. The capacity for refreezing of meltwater within the snowpack is simulated based on the temperature and mechanical properties (e.g. density, pore volume) (e.g. Illangasekare and others, 1990; Mölg and others, 2008). Changes in the density and temperature due to refreezing are distributed evenly through the snowpack and are determined by a weighted average. Storage capacity and densification of the snowpack are simulated using the approach of Vionnet and others (2012). Further, settling and compaction of snow is simulated using the viscosity approach (Sturm and Holmgren, 1998). The key model parameters and approaches are presented in Table 1.

Table 1. Input parameters used in the fully distributed surface energy and mass-balance model

3.2. Model evaluation

The mass-balance model is evaluated using distributed simulations over the 10-year meteorological record between June 2010 and November 2020. To capture the entire glaciological mass balance in the first year (2010/11) meteorological data for the period 29 March–25 June 2010 are reconstructed using meteorological data over the same period from the following mass-balance year (2011/12). The parameterisations applied in the model (e.g. precipitation scaling, cloud metrics, air temperature lapse rate, surface roughness lengths, albedo for snow, ice and firn, density of snow and ice, snow/rain threshold, subsurface routines) have been used in a large number of energy and mass-balance studies (e.g. Mölg and others, 2008, 2009a, 2009b; Gillett and Cullen, 2011; Mölg and Kaser, 2011; Conway and Cullen, 2013, 2016; Conway and others, 2015; Cullen and Conway, 2015). Of particular relevance, Conway and Cullen (2013) determined roughness lengths of momentum, temperature and humidity for Brewster Glacier and a Monte Carlo analysis was used to evaluate the uncertainties in calculating the turbulent heat fluxes using the bulk aerodynamic method. Conway and others (2015) developed cloud metrics and parameterisations to quantify the influence of clouds on incoming shortwave and longwave radiation, while the albedo parameterisation of Oerlemans and Knap (1998) is optimised by Conway and Cullen (2016). These findings were used to guide the current approaches to optimise model parameterisations.

We conducted a series of runs to establish the optimal parameter set and to quantify some of the uncertainties due to parameter choice. A comprehensive analysis enabled nine model runs to be selected (Table S1). Initial runs show that the precipitation correction (Section 3.1) is overestimated. The chosen model configurations are based on the expected influence of precipitation, lapse rate, rain/snow threshold, roughness length and albedo on mass balance on Brewster Glacier. To analyse model performance, we compare the glaciological observations of glacier-wide winter (Fig. S1), summer (Fig. S2), annual (Fig. S3) and altitudinal annual (Fig. S4) mass balances to modelled output. For each of the nine runs (A–I), we calculate the correlation coefficients (R), root mean square error (RMSE) and mean bias error (MBE). The model runs are ranked (1 – highest R, lowest RMSE and MBE and 9 – lowest R, highest RMSE and MBE). The three model runs which ranked the lowest are selected and further evaluated against the winter, summer and annual calibrated glaciological observations (Fig. S5). From this analysis, MOD _I (hereby referred to as MOD _MB) is selected as the baseline run in this study.

3.3. Nomenclature for comparing glaciological and modelled mass balance

To clarify the terms used to compare glaciological and modelled mass balances, a schematic of the annual cycles of hypothetical mass balance for the terminus, mid- and high-elevation sections of Brewster Glacier is presented in Figure 2. Periods of winter and summer mass balance captured by the glaciological observations are referred to as Obs _MB, while modelled mass balances over the same glaciological period are defined as MOD _{observed-like dates}. The average dates of the summer and winter glaciological mass balance are shown as vertical lines in Figure 2, but it should be noted these dates vary each year since a floating-date system is used (Cullen and others, 2017). The modelled seasonal mass balance that explicitly accounts for the duration and extent of glacier-wide ablation and accumulation is defined as MOD _{full period}. We refer to the difference between the glaciological observational period (Obs _MB and MOD _{observed-like dates}) and MOD _{full period} as the ‘missing mass balance’ (shown in yellow, Fig. 2). We refer to the calibrated (or reconstructed) glaciological mass-balance record of Brewster Glacier developed by Sirguey and others (2019) as Obs _{MB,calibrated}, which is compared to MOD _{full period} to quantitatively assess modelling uncertainties in mass balance.

Fig. 2. Schematic highlighting the glaciological mass balance (Obs _MB), periods of missing mass balance and modelled periods at the terminus, mid- and high-elevation regions of Brewster Glacier. The approximate period when end of summer and winter glaciological observations are obtained is shown, as are the range of the variable start dates. The full definitions of the acronyms used to define the different model simulations are provided in the text.

4. Results

4.1. Modelled glacier mass balance

The modelled and observed glacier mass balance during the glaciological period, and the performance of MOD _MB at seasonal and annual glacier-wide and annual altitudinal scales are shown in Table 2.

Table 2. Modelled and observed glacier mass balance

Model performance using the correlation coefficient (R), root mean square error (RMSE) and mean bias error (MBE) at glacier-wide and altitudinal scale is also shown. Mass balances are shown for the glaciological period (2010/11 to 2019/20).

The winter mass balance derived by the mass-balance model slightly underestimates the winter glaciological mass balance. On average, the results indicate that the mass-balance model reproduces glaciological mass balance in winter quite well, while modelled ablation in summer is more negative for the observed-like-dates (Table 2). The altitudinal variations in observed and modelled mass balance demonstrate that the winter balances at the higher elevations are not as positive in the model run as compared to observations, which is compensated for by a more positive winter mass balance at the lower elevations (Fig. 3). The failure to reproduce annual balance in the lowest part of the ablation area is mainly controlled by an overestimate in the winter mass balance that persists in to summer. Yearly variations in the altitudinal mass balances are observed for both Obs _MB and MOD _{observed-like dates}, but modelled variations are somewhat less pronounced than observed, particularly for summer mass balance, though we note that glaciological mass balance >2000 m is extrapolated using observations from below this elevation threshold, so are subject to a larger uncertainty.

Fig. 3. Altitudinal mass balance of Brewster Glacier over the study period. (a) Annual balance; (b) winter balance; (c) summer balance; and (d) hypsometry. The red and grey shading depicts one std dev. from the observed glaciological mass balance (Obs _MB) mean and the modelled mass-balance (MOD _{observed-like dates}) mean, respectively.

4.2. The missing mass balance

Figure 4 enables the annual cycles of modelled ablation and accumulation over the 10-year period to be assessed and compared to the end and start dates of the glaciological mass-balance floating-date system. While differences are observed in the modelled end of ablation season for each of the nine evaluated model runs (Fig. S6), it is important here that we emphasise that this pattern of an extended ablation season is consistent across each run. Figure 4 indicates that the ablation season (summer balance) typically extends beyond the glaciological end of summer observation window (12 March to 8 April). Over the 10 years, the end of the modelled ablation period, based on minimum mass balance, ranges between 8 April and 13 May. This indicates that the end of summer glaciological observations are consistently missing melt and are largely responsible for the known positive bias in the glaciological mass balance (Sirguey and others, 2019). The modelled end of the accumulation period (winter balance) ranges between 4 November and 29 November, while the observation window is 30 October to 1 December. The mean modelled glacier-wide ablation season (summer balance) is 166 days on Brewster Glacier, while the accumulation season (winter balance) is longer (199 days).

Fig. 4. (a) The modelled mass balance (MB) for each year over the study period (a mass balance year is presented as 1 April–31 March); (b) snow accumulation (M _acc); and (c) surface melt (M _melt) during each of the mass-balance years. The red line represents the mean mass balance over the 10 years. The red vertical lines indicate the range of dates for the end of the ablation period (summer balance) and the end of the accumulation period (winter balance) for each of the modelled mass-balance years.

To explore this further, the difference between the observed glaciological mass balance date at the end of summer and the end of the modelled ablation period (end of modelled summer balance) is analysed (Fig. 5a). The timing of when the modelled accumulation season (winter balance) begins for different elevations across the glacier for each year is presented in Figure 5b. It is clear that mass loss is being missed at all elevations across the glacier due to the glaciological observations being obtained prior to the end of the ablation season, with the largest mass loss variability in the lower ablation zone. As expected, the transition from ablation to accumulation varies with altitude (Fig. 5b), with a degree of inter-annual variability. For example, large amounts of missing mass balance are seen in 2011, 2015 and 2017, while 2014 and 2018 are observed to have less missing mass loss. On average, the mass loss missed by the glaciological observations is equal to −516 mm w.e. a⁻¹, which represents 35% of the modelled annual mass balance over the 10-year period. Figure 6 demonstrates further that the largest proportion of the missing mass balance occurs at the lower elevations of the glacier (below 2000 m).

Fig. 5. (a) Mass balance of each year for the period between the observed glaciological mass-balance date at the end of summer and the end of the modelled ablation period (end of modelled summer balance). The mean mass balance for this period of difference between observed and modelled summer balance is shown in red; and (b) the date of the start of the modelled accumulation season (winter balance) for each year at different elevations of the glacier.

Fig. 6. Average spatial differences between modelled annual balance using the glaciological observation periods (MOD _{observed-like dates}) and modelled mass balance (MOD _{full period}) not constrained by the floating-date system. The units are mm w.e., with the largest differences observed in the ablation zone (bottom of figure).

A comparison of the observed glaciological mass balance and calibrated glaciological mass balance (Obs _{MB,calibrated}; Sirguey and others, 2019) records to the modelled mass balance demonstrates clearly that the model generally captures the year to year variability well (Fig. 7). On an annual scale, an RMSE of 598 mm w.e. and MBE of 146 mm w.e. are observed between the modelled mass balance and the calibrated glaciological mass balance. The winter balance demonstrates an RMSE of 304 mm w.e. and MBE of 267 mm w.e., while an RMSE of 587 mm w.e. and an MBE of −121 mm w.e. are observed between the modelled and calibrated glaciological summer balance. It is evident that both the modelled and calibrated annual balances are negative each year (Fig. 7b), which is not the case for the glaciological and modelled annual balances using the floating-date system (Fig. 7a). The evidence from the mass-balance modelling supports the findings of Sirguey and others (2019), which is that Brewster Glacier has not observed any positive mass balances in the period 2010–2020 and is in a state of significant decline.

Fig. 7. The glacier-wide seasonal and annual mass balance (MB) of (a) Obs _MB and MOD _{observed-like dates}; and (b) Obs _{MB, calibrated} and MOD _{full period}.

4.3. Winter melt

Figure 8 shows the interannual variability in the total amount of ‘missing mass balance’ as well as the amount of modelled melt that occurs during winter. Interannual variability in melt is observed during both periods. The winter of 2014 experiences the highest amount of surface melt, with almost 23% of the snow accumulation lost to surface melting across the glacier. On average, across the glacier, 2–32% of the snow accumulation during the accumulation season melts (Fig. 9a).

Fig. 8. Total surface melt that occurs in the period between the end of the summer glaciological mass-balance observations (end of observed ablation season) and the end of the modelled ablation period, as well as the amount of modelled melt during the accumulation season (winter period).

Fig. 9. Spatial patterns of (a) percentage of snow accumulation lost through surface melt during the modelled accumulation season (winter); and percentage frequency of melt days during (b) June–August; and (c) December–February. Melt days are assumed when daily average Q _M > 0.

At low elevations, almost 19% of the days experience melting (defined as days where Q _M > 0) during June–August, while this reduces to 2% at higher elevations (Fig. 9b). On the other hand, during December–February, this value is notably higher with 99 and 94% at lower and higher elevations, respectively (Fig. 9c). These values represent the averages over the 10 years and both spatial and temporal variations are observed.

4.4. Vertical gradients of mass fluxes

The spatial and temporal variations of the mass fluxes indicate that higher melt and lower accumulation led to the high mass loss observed during 2017/18 (Fig. 10). Similarly, reduced accumulation and refrozen meltwater, and higher than average subsurface melt and sublimation contributed to a negative annual mass balance in 2011/12. Refrozen meltwater is an important mass flux at Brewster Glacier. On average, ~10% of surface and subsurface melt is refrozen over the study period. While 434 mm w.e. of surface and subsurface meltwater is refrozen on an annual basis, this value increases with elevation (average difference of 427 mm w.e. between the lowest and highest elevation points). This suggests that meltwater and possibly rain is often being refrozen within the snowpack. Overall, the contributions of surface deposition and sublimation are of similar magnitude, however larger variations with elevation are observed for surface deposition.

Fig. 10. Altitudinal profiles of mass-balance components for each year over the study period. (a) Annual mass balance (MB); (b) snow accumulation (M _acc); (c) refrozen meltwater (M _ref); (d) surface deposition (M _dep); (e) surface melt (M _melt); (f) surface sublimation (M _sub); and (g) subsurface melt (M _subsurfmelt). The mean annual profiles are shown in red. The components of each mass-balance year run from 1 April to 31 March. The negative mass fluxes represent the fluxes that contribute to mass loss while the positive fluxes are associated with mass gain.

4.5. Surface energy balance

Net radiation (Q _Net) dominates the surface energy balance during the study period. Given average net longwave radiation only varies moderately each month, net radiation is governed by the distinct seasonal variability in net shortwave radiation (Tables 3 and 4). The average annual net shortwave radiation is 38 W m⁻², while the average net energy loss by net longwave radiation is −25 W m⁻². Annually, sensible heat flux is of equal magnitude to net longwave radiation and provides 25 W m⁻² of energy to the glacier. While the latent heat flux only supplies an average of 2 W m⁻² of energy to the surface annually, it becomes an important energy source for ablation during summer. Heat from rain is the smallest energy source to the surface, contributing as little as 1 W m⁻² annually, but supplies more energy during summer during heavy rainfall events. The negative penetrating shortwave radiation creates subsurface melt during summer as the conductive heat flux is small. The subsurface melt is small during winter when the conductive heat flux increases. Annually, penetrating shortwave radiation removes −6 W m⁻² of energy from the surface to the subsurface, while the conductive heat flux contributes 8 W m⁻² of energy to the surface energy balance. The main melt period extends from November–April, with an average of 80 W m⁻² during this period. Annually, melt energy is 43 W m⁻², with marked increases during periods of melting driven by increases in net shortwave radiation, sensible heat, latent heat and heat from rain (annual average when Q _M > 0 is 117 W m⁻²) (Table 4).

Table 3. Mean monthly values of the energy fluxes over the study period (2010/11 to 2019/20)

Values are in W m⁻².

Table 4. Annual average of the energy fluxes over the study period (2010/11 to 2019/20)

Square brackets give values during periods of surface melt (Q _M > 0). Values are in W m⁻².

Subsurface melt is not included in Q _M.

Altitudinal profiles of the surface energy fluxes within 100 m elevation bands during all periods and periods of melting (Q _M > 0) are shown in Figure 11. The percentage of periods when melting occurs across the glacier (Q _M > 0) is also shown (Fig. 11c). Net shortwave radiation governs variability in net radiation, which on a glacier-wide scale contributes to 59–65% of the energy available for melt (Q _M) (Table 4). The turbulent heat fluxes are also important contributors to the energy balance, with sensible heat flux responsible for 32–33% of melt energy and latent heat flux, 10–15%. The significance of the conductive heat flux decreases during periods of melt across the glacier, while heat from rain varies only slightly, both spatially and during both periods (average is 3–4% for both periods). The penetration of shortwave radiation to the subsurface of the glacier is important. During periods of melt, penetrating shortwave radiation leads to 7–15% of the energy surplus at the surface being lost to the subsurface layer and is most important at the lower elevations.

Fig. 11. Altitudinal profiles of average energy fluxes. (a) All periods over the study period; (b) periods of melt when Q _M > 0; and (c) percentage of periods of melt when Q _M > 0 for each mass-balance year. Fluxes are calculated at 100 m elevation bands.

4.6. Surface energetics during the main melting period

To better understand the controls on the variability of melt energy, the correlation between the energy fluxes, melt energy and air temperature during the main melting period (assumed to November–April) is examined (Table 5). Even though net shortwave radiation contributes a large proportion of the energy for melt it has a weak and not statistically significant (p < 0.05) relationship to melt energy (R = 0.49). Instead, the variability of melt energy is mainly controlled by the magnitude of the sensible heat flux (R = 0.88), latent heat flux (R = 0.68), heat from rain (R = 0.69) and net longwave radiation (R = 0.66), highlighting the importance of airmass variability on glacier melt in the Southern Alps. Air temperature exhibits a strong and significant correlation to melt energy (R = 0.83). Air temperature is also significantly correlated to the sensible heat flux (R = 0.71) (air temperature is used to estimate the sensible heat flux) and net shortwave radiation (R = 0.70), while weaker correlations that are not statistically significant are observed for latent heat flux (0.42), heat from rain (0.39) and net longwave radiation (R = 0.27).

Table 5. Correlation coefficients (R) of the main ablation seasonal anomalies (November–April) of the energy fluxes and melt energy (Q _M) and air temperature (T _a) over the study period (2010/11 to 2019/20)

Statistical significance is represented in bold (p < 0.05).

4.7. Influences on the interannual variations in mass balance

The statistical relationship between snow accumulation and the energy available for melt and seasonal anomalies of moisture-related atmospheric variables (wind speed, relative humidity, precipitation and cloud cover) are shown in Table 6. The wind data at AWS_lake are influenced by both local and regional scale atmospheric processes, with the local winds in particular modified by the surrounding topography. Previous studies have shown that despite the influence of cloud cover on the magnitude of the individual energy fluxes, average daily values of melt energy are similar under both clear-sky and overcast conditions, mainly due to the ability of melt to be maintained over longer periods under overcast conditions (Conway and Cullen, 2016). Precipitation shares a strong correlation with snow accumulation (R = 0.67), while wind speed and cloud cover are moderately correlated to snow accumulation. The evidence of a statistical relationship between mass gain and energy available for melt and moisture-related atmospheric variables suggests that variability in mass balance is in part influenced by air mass variability.

Table 6. Correlation coefficients (R) of the main ablation seasonal anomalies (November–April) of moisture-related atmospheric variables (wind speed (U), relative humidity (RH), precipitation (P) and cloud cover) and melt energy (Q _M) and snow accumulation (M _acc) during the study period (2010/11 to 2019/20)

Statistical significance is represented in bold (p < 0.05).

Modelled glacier mass balance indicates that interannual variability in winter mass balance is primarily influenced by changes in snow accumulation (R = 0.98), with contributions from refrozen meltwater (R = 0.61) (Table 7). Summer mass balance is driven primarily by changes in surface melt (R = −0.98) and subsurface melt (R = −0.72), with contributions from snow accumulation (R = 0.86), refrozen meltwater (R = 0.76) and surface deposition (R = −0.74).

Table 7. Correlation coefficients (R) of the seasonal anomalies (both accumulation and ablation) of the mass fluxes and glacier mass balance (MB)

Statistical significance is represented in bold (p < 0.05).

Snow accumulation and liquid precipitation (rain) during modelled accumulation and ablation periods are presented in Table 8. Interannual variability in snowfall and rainfall is observed for winter and summer. On average, almost 71% of the total precipitation falls as snow during winter. During summer, 12% of the total precipitation is snowfall, with the highest level of summer accumulation observed during 2019/20. During 2019/20, 430 mm w.e. of snowfall is observed, which was likely responsible for the extended period of accumulation observed in Figure 4a. Further, some of the mass-balance years with higher rainfall and reduced snow accumulation during summer (e.g. 2015/16 and 2018/19) exhibit more negative summer mass balances (Table 8 and Fig. 7b). The relationship between summer snowfall and summer mass balance exhibits a strong and significant correlation (R = 0.86), indicating that like other mid-latitude glaciers (e.g. Oerlemans and Klok, 2004), summer accumulation can play an important role in interannual mass-balance variability on Brewster Glacier. Finally, Table 8 clearly shows that proportionally, more rain falls in winter than snow falls in summer, which is perhaps the starkest evidence that both winter and summer are becoming less effective at storing mass.

Table 8. Snow accumulation and liquid precipitation (rain) during modelled accumulation and ablation season

5. Discussion

5.1. Reconciliation of glaciological mass balance

Glaciological in situ mass-balance measurements are strongly influenced by the time periods associated with the mass-balance quantities and the spatial extent of the sampling scheme (Mayo and others, 1972; Jansson, 1999). While the glaciological measurement campaign on Brewster Glacier consists of extensive geo-located probing of snow depth across the glacier and ablation stakes drilled primarily along the central flowline, the spatial density of the sampling network and the difficulty in measuring accumulation and ablation in the upper accumulation areas (Fig. 1) may explain some of the discrepancy between the glaciological and geodetic assessment of mass loss. The altitudinal profiles of modelled mass balance suggest that the glaciological winter balance overestimates accumulation above 2000 m a.s.l., a region of the glacier that is under-sampled (Cullen and others, 2017). It is important here that we also highlight some of the uncertainties in the processes simulated by the model that may account for some of the differences observed between the seasonal modelled and glaciological observations (Figs 3 and 7). For example, the model overestimates the winter balance in the lower elevation regions (average difference of 967 mm w.e.), which persists through to summer resulting in overestimation of the modelled annual mass balance at low elevations (Fig. 3). It is possible that some of these differences are associated with the input parameters used in the model (Table 1). While the input parameters are based primarily on observations and calculations that are specific to Brewster Glacier, some parameters including the rain/snow threshold and penetrating shortwave radiation scheme have not been validated for Brewster Glacier and are subject to uncertainties.

Modelled glacier mass balance shows an average mass loss of −516 mm w.e. a⁻¹ is not captured by glaciological observations. While mass loss is observed across all elevations, this missing mass balance is primarily associated with additional ablation at lower elevations. The modelled glacier mass balance suggests that Brewster Glacier is experiencing an extended ablation season, which may have implications on the end of summer snowline survey, which is conducted in mid-March. It is clear that melting extends beyond this period at all elevations of Brewster Glacier. Glaciers, globally, are experiencing dramatic changes in response to a warming climate, with the extension of the ablation season observed across several regions (Thibert and others, 2018; Marshall, 2021). This extended ablation season results in the expansion of the ablation region and increased melt due to longer duration of the lower albedo ice-exposed surface (Marshall, 2021). While in some years the modelled end of the ablation season is almost uniform across Brewster Glacier, for other years large variations in timing are observed across different elevations (Fig. 5b). This finding highlights the complexity in determining the end of the ablation season, particularly for the end of summer snowline aerial survey, which must be done before snowfall commences in the higher elevation sections of the glacier. Furthermore, the extended ablation season raises key questions around mass-balance reconstruction from end of summer snowline surveys. Despite the timing issues, strong significant correlations are observed between the annual modelled mass-balance and glaciological observations (R = 0.74), and the EOSS elevations (R = −0.90). Similarly, the winter and summer glaciological mass balances also indicate significant correlations with modelled winter (R = 0.86) and summer mass balance (R = 0.74). This suggests that both the glaciological observations and EOSS provide very useful information in determining glacier mass balance in the Southern Alps, but should be calibrated with other techniques including geodetic estimates and numerical models to account for any biases and uncertainties.

The discrepancies between the glaciological observations and geodetic data suggest that a homogenised mass-balance record is required for Brewster Glacier. While the glaciological method is generally assumed to be the most reliable and an important baseline for monitoring glaciers, this method is also subject to uncertainties, biases and systematic errors. Previous studies have demonstrated the importance of calibrating and combining glaciological observations with geodetic volume changes and distributed numerical models (e.g. Machguth and others, 2006; Cogley, 2009; Huss and others, 2009; Marinsek and Ermolin, 2015; Thomson and others, 2017), which produces better temporal and spatial representation of mass-balance changes and a more homogenised time series of mass-balance records. While distributed mass-balance modelling is an important tool for long-term mass-balance monitoring (Machguth and others, 2006), this is yet to be implemented for Brewster Glacier on a long-term basis. The optimisation and evaluation of the distributed mass-balance model with observational data in this study is the first step towards the full homogenisation of the mass balance for Brewster Glacier (Huss and others, 2009). It has not been attempted in this study as we did not want to introduce additional input uncertainty into the present modelling framework.

5.2. Surface energetics and drivers of melt: implications and comparison with previous studies

The partitioning of the surface energy balance of Brewster Glacier is similar to other mid-latitude glaciers (Greuell and Smeets, 2001; Fitzpatrick and others, 2017), with net radiation driven primarily by net shortwave radiation dominating the available melt energy (63%). The turbulent sensible and latent heat fluxes are important contributors to melt providing the majority of the remaining energy to the surface (32 and 14%, respectively). This highlights the need to resolve any uncertainties in the turbulent heat flux parameterisation scheme, including the spatio-temporal variations in the surface roughness lengths (Conway and Cullen, 2013; Fitzpatrick and others, 2017). The heat from rain plays a minor role, contributing 3% to the surface energy balance. Penetrating shortwave radiation removes 12% of the surface energy to the sub-surface layer. The dominant control of net radiation on the surface energy balance and strong albedo effect is observed during 2016/17 and 2017/18, with the difference in melt energy largely driven by changes in net radiation (Table 4).

While net radiation is a key component of the surface energy balance, variability in melt energy is strongly controlled by changes in the turbulent heat fluxes, heat from rain and net longwave radiation (Table 6). Despite the relatively small contribution of heat from rain to the surface energy balance, its significant correlation to melt energy suggests heat from rain is potentially an important energy flux during extreme rainfall events, which are projected to become more frequent in a warming climate (Intergovernmental Panel on Climate Change, 2021; Bodeker and others, 2022). A significantly strong correlation is observed between melt energy and air temperature (R = 0.83). However, variability in melt energy is not only associated with air temperature changes, with atmospheric moisture reflected through precipitation, cloud cover and wind speed accounting for some of the observed variability in melt energy and snow accumulation (Table 6). As demonstrated in previous studies (Gillett and Cullen, 2011; Conway and Cullen, 2016; Abrahim and others, 2022), melt energy and glacier mass-balance variability in the Southern Alps cannot be solely explained by changes in air temperature, with changes to airmasses and low-level moisture advection playing an important role in this region.

Past studies have shown that the interannual variability in the mass balance is primarily driven by summer mass balance (Anderson and others, 2010; Cullen and others, 2017). To explore this further, the processes that contribute to both summer and winter mass balance are examined, with a particular focus on snowfall variability and refrozen meltwater. On average, at low elevations, 32% of the snow accumulation during winter is lost through surface melt, while at higher elevations only 2% is removed. Furthermore, interannual variability in winter snow accumulation is observed, with the two most positive winter mass years (2012/13 and 2018/19) experiencing the largest quantity of winter snowfall (Table 8). In the Southern Alps, large snowfall events have been shown to contribute as much as 35–40% of the total annual snow accumulation (Porhemmat and others, 2021). The rain/snow threshold and the magnitude of extreme snowfall events will therefore have implications on glacier mass balance in the Southern Alps as the climate continues to warm and extreme events become more frequent in the future. Like snowfall variability, refrozen meltwater is an important mass flux that influences both winter and summer mass balances. On average, at least 10% of the melted water on Brewster Glacier is being refrozen. While the glaciological mass balance provides an estimation of refreezing in the annual surface layer, the modelled quantity of refrozen meltwater includes refreezing of percolating meltwater through the snow and firn layers. Internal accumulation in firn due to refreezing can lead to an underestimation of the annual mass balance of a glacier and is therefore a potential source of systematic error if not included in mass-balance estimations (Schneider and Jansson, 2004). Refrozen meltwater appears to act as a net buffer to mass loss across all elevations on Brewster Glacier and rapid retreat could occur if this refrozen meltwater is lost, rather than being stored in the system, as retreat is exacerbated in the future.

Although no significant correlation is established between the modelled winter and summer mass balance at Brewster Glacier (R = −0.14; p = 0.70), we find that summer mass balance is significantly correlated to the annual glacier mass balance (R = 0.89; p = 0.0005). This highlights the importance of melt processes during summer on the overall annual mass balance of Brewster Glacier. A strong, significant correlation between modelled summer snow accumulation and summer mass balance is also observed (R = 0.86). This suggests that snowfall during summer reduces melt due to the higher albedo of fresh snow, which results in less absorbed solar radiation (Oerlemans and Klok, 2004). This reduced surface melt leads to greater variability in mass balance during summer (Kuhn, 1993), which is observed in this study. Some of the less negative summer balance years (e.g. 2016/17) were characterised by a larger proportion of summer snowfall (Table 8). It has been previously suggested by Oerlemans and Klok (2004) that two large snowfall events could compensate for a 1 K higher mean temperature during the ablation season at a mid-latitude glacier site. The findings of this study are consistent with findings on other mid-latitude glaciers (Marshall and Miller, 2020), demonstrating that summer snowfall plays a key role in protecting glaciers in summer. Thus, changes in the amount of snowfall during the ablation season will likely play an important role in the future fate of glaciers in the Southern Alps.

5.3. Sensitivity and uncertainties in the model

The sensitivity of Brewster Glacier to changes in air temperature and precipitation is examined in Figure 12 to enable a comparison to other studies exploring glacier response to a simple climate forcing (e.g. Anderson and others, 2010). This modelling exercise demonstrates that the response of mass balance to changes in air temperature is not symmetrical, with a +1 K in temperature resulting in a lowering of the mean glacier-wide mass balance by −1639 mm w.e. a⁻¹, while a -1 K decrease in air temperature results in a mass gain of 1451 mm w.e. a⁻¹. This sensitivity to changes in air temperature increases dramatically at lower elevations on the glacier (Fig. 12). An increase from +1 to +2 K results in an average glacier-wide mass-balance difference of −1883 mm w.e., while the difference between the −1 and −2 K simulation is 1230 mm w.e. The large sensitivity of mass balance to changes in air temperature is likely associated with variations in surface albedo, partitioning of rain/snow and the length of the ablation season (Oerlemans and Hoogendoorn, 1989). An increase in precipitation of 10% leads to a mass gain of 477 mm w.e. a⁻¹, while a 10% precipitation decrease results in a mass loss of −476 mm w.e. a⁻¹. The largest effect of increases in precipitation, seen more clearly in a 20% increase, is observed at the higher elevations, where the proportion of solid precipitation and its associated albedo feedback increases. The moderate changes in mass balance for a 10% precipitation change suggest that large increases in precipitation would be required to offset the mass loss associated with the expected atmospheric temperature rise over the 21st century. The estimates of temperature sensitivity here are lower than those reported by Anderson and others (2010) (−2000 mm w.e. °C⁻¹ a⁻¹), while precipitation sensitivity here is almost 20% higher than reported by Anderson and others (2010) (400 mm w.e. a⁻¹ for 10% precipitation change). The findings of this study confirm that Brewster Glacier, like other mid-latitude glaciers, remains very sensitive to changes in air temperature.

Fig. 12. Sensitivity of the annual mass balance to (a) air temperature; and (b) precipitation.

Redistribution of snow by wind and avalanching can result in modified snow accumulation on glacier surfaces (Mott and others, 2008; Dadic and others, 2010; Sauter and others, 2013). Several studies have shown the key role snow redistribution and preferential deposition by wind play in mountainous glaciated regions in the European Alps (Mott and others, 2008; Dadic and others, 2010). The distributed mass-balance model used in this study did not consider snow redistribution, which may have implications on modelled mass balance. Interactions with the local topography and wind dynamics can lead to spatial variations in snow depths, with exposed windward areas more likely to exhibit less deposition (Liston and Sturm, 1998). The lower snow accumulation observed below 1800 m a.s.l. in the glaciological winter balance (Fig. 3b) suggests that local-scale processes beyond those imposed by hypsometry and not simulated by the model may be important in explaining glacier mass-balance variability. Snow deposition from avalanching at ~2000 m a.s.l. associated with steep terrain, which is observed in the glaciological winter balance (Fig. 3b), is also not captured by the distributed mass-balance model. Furthermore, elevational gradients are applied to distribute meteorological variables (e.g. temperature lapse rate) across the glacier. The AWS data are collected at 1650 m a.s.l. and are obtained over bedrock, which may have implications on the daily cycles of air temperature that are not fully accounted for (Cullen and Conway, 2015). Additionally, a constant lapse rate is used to spatially distribute air temperature, which leads to larger uncertainties with elevation. Using a constant lapse rate for day and night may fail to capture any surface cooling during the night and could possibly lead to overestimation of melt. Unlike air temperature, no gradients for wind speed, precipitation or cloud cover are applied and any uncertainties in these variables may extend across the entire glacier. The calculation of the turbulent heat fluxes is based on both observational findings and assumptions. While the surface roughness lengths used are those observed on Brewster Glacier using eddy covariance measurements, surface roughness lengths are known to exhibit both spatial and temporal variability (Conway and Cullen, 2013). Variations in temperature and vapour pressure gradients can affect the spatial and temporal variability of surface roughness lengths and the roughness values derived in Conway and Cullen (2013) may not be representative over extended time periods or across the entire glacier.

Intensity and phase of precipitation at Brewster Glacier also has limitations as precipitation is scaled from a low elevation site (Cullen and Conway, 2015). The results presented here suggest the precipitation record deserves to be further scrutinised. The rain/snow threshold for Brewster Glacier also remains uncertain and new methods to better understand the physical processes governing precipitation in the Southern Alps are still required. The rain/snow threshold is a key factor in glacier mass balance, influencing both snow depth and snow water equivalent (Loth and others, 1993) and has been shown to exhibit both spatial and temporal variations in the Northern Hemisphere (Kienzle, 2008; Jennings and others, 2018). Snow accumulation is an important mass flux at Brewster Glacier, so the magnitude and phase of precipitation are critical in controlling variability in glacier mass balance.

6. Conclusion

A fully distributed mass-balance model has been used to assess seasonal variability in the surface energy and mass balance of Brewster Glacier in the Southern Alps of New Zealand over a 10-year period. The approach taken was to compare output from the mass-balance modelling to a time series of mass balance revised using geodetic methods to assess uncertainties in the glaciological mass balance. The approach has also enabled the atmospheric controls on seasonal mass balance to be carefully assessed and the sensitivity of Brewster Glacier to our present climate to be tested. The main findings and conclusions of this study are:

1. (1) The observed glaciological summer balance has been overestimated on Brewster Glacier due to the failure to capture mass loss during the latter part of the ablation season. Mass-balance modelling revealed that, on average, 38 days of ablation have been missed by the mass-balance monitoring programme. The average modelled length of the ablation season over Brewster Glacier in the present climate is 16 November to 30 April (166 days), with the accumulation season accounting for the remainder of the year (1 May to 15 November; 199 days).

2. (2) The average annual mass loss not accounted for by the glaciological observations is −516 mm w.e., which represents 35% of the modelled annual mass balance. Annual glacier-wide mass balance for Brewster Glacier is estimated to be −1481 mm w.e. a⁻¹ over the period 2010–2020, with no individual year experiencing a positive annual balance.

3. (3) Mass-balance modelling shows that summer balance is significantly correlated to annual balance, however, ablation processes are not only important in summer. Even at high elevations, melting during the accumulation season is not uncommon on Brewster Glacier, with 2–32% of the snowfall melting during the accumulation season. As warming continues, rain-on-snow events are likely to become more common in winter and the observed loss of seasonal snow will likely further compromise winter balance on Brewster Glacier.

4. (4) Refrozen meltwater is an important mass accumulation flux at Brewster Glacier, with 10% of surface and subsurface melt being refrozen. The role frozen meltwater plays in mass balance of other glaciers in the Southern Alps remains largely unknown but evidence from this study is that it warrants further investigation.

5. (5) Net radiation, driven primarily by net shortwave radiation, is the dominant source of energy, contributing to 63% of the total energy available for melt on a glacier-wide scale. The penetration of shortwave radiation into the subsurface of the glacier during periods of melt is important. This energy is lost to the subsurface (−12%) and does not contribute directly to surface melt. The turbulent sensible (32%) and latent heat (14%) fluxes account for most of the remaining melt energy (46%). On average, the energy associated with precipitation is a small contributor to melt (3%) but becomes proportionally larger during rain-on-snow events.

6. (6) There is a strong correlation between interannual variations in melt energy and air temperature, which supports findings that Brewster Glacier, like other mid-latitude glaciers, are most sensitive to perturbations in temperature. However, the variability in melt energy during the ablation season, in particular the melt extremes, is controlled by an increase in energy from the turbulent sensible and latent heat fluxes, net longwave radiation and the heat flux from precipitation.

7. (7) A strong correlation between snow accumulation in summer and summer balance indicates that solid precipitation events during the ablation season do reduce melt effectively by adding mass and lowering the amount of absorbed solar radiation. The sensitivity of summer snowfall on Brewster Glacier mass balance has not been fully tested but the observed correlation between snowfall and summer balance could be explored further by determining the sensitivity of glaciers in the Southern Alps to extreme precipitation events, particularly as they are projected to become more frequent in the future.

Overall, this study provides new insights about the atmospheric controls on mass balance on Brewster Glacier, both spatially and temporally in a way that has not been captured before. We show that the glacier–climate relationship is complex and even though air temperature and precipitation are the main controls on glacier mass balance, atmospheric moisture fluxes play a key role in the seasonal variability of ablation and accumulation. The mass balance on Brewster Glacier is clearly negative at present, with the evidence presented here supporting other research efforts that mass loss during the observational period of glaciers in the Southern Alps has never been greater than what is occurring presently. The glaciological mass balance requires a complete re-assessment using an approach similar to Huss and others (2009) given the amount of mass loss that is not being captured at the end of the ablation season. It has not been attempted in this study due to our reluctance to introduce additional input uncertainty into the present mass-balance modelling framework. The uncertainty in the glaciological mass balance also raises questions about how to most effectively use the EOSS aerial observations as a proxy for glacier mass balance more widely across the Southern Alps (Lorrey and others, 2022). Repeated regular geodetic estimates are becoming more readily available, providing a robust approach to compensate for the propagation of biases found in other approaches to estimate mass balance. The challenge that lies ahead is to build a flexible and multi-proxy approach to accurately document the unprecedented changes the unique glaciers in the Southern Alps are currently undergoing.