Introduction
Information about the ‘climate sensitivity’ of glaciers is useful for predicting the terminus response, meltwater runoff and glacier contribution to sea-level rise under future climate change (Reference OerlemansOerlemans and others, 1998), and for using palaeoglacial extents to reconstruct past climate changes (Reference Anderson and MackintoshAnderson and Mackintosh, 2006). The climate sensitivity, or mass-balance change for a given temperature or precipitation change, has been calculated using simple empirical schemes (Reference Oerlemans and FortuinOerlemans and Fortuin, 1992; Reference Braithwaite, Zhang and RaperBraithwaite and others, 2002) founded on measurement and modelling of energy and mass balance of glaciers (Reference OerlemansOerlemans, 1992). An important finding of these studies is that climate sensitivity varies geographically. In particular, high-precipitation glaciers are more sensitive than glaciers in more arid regions.
The glaciers of the Southern Alps of New Zealand have received less attention than other areas, such as the European Alps. However, these glaciers are important because of their Southern Hemisphere maritime setting, and the high precipitation they experience, resulting in glaciers that are likely to be highly sensitive to climatic change. Although the complete loss of the glaciers of New Zealand would result in a negligible (0.15 mm) contribution to eustatic sea-level rise, they provide insight into the mass-balance response of maritime glaciers in general, which are predicted to make a major short-term contribution to sea-level rise (Reference SolomonSolomon and others, 2007), through both mass-balance changes and dynamic response.
Of the glaciers for which mass-balance sensitivity has been calculated, Franz Josef Glacier has one of the highest sensitivities (Reference Oerlemans and KnapOerlemans and Knap, 1998; Reference Braithwaite, Zhang and RaperBraithwaite and others, 2002; Reference De Woul and HockDe Woul and Hock, 2005; Reference Anderson, Lawson, Owens and GoodsellAnderson and others, 2006). This is because the western flank of the Southern Alps experiences a hyper-maritime climate, with annual precipitation >12 m in places (Reference Griffiths and McSaveneyGriffiths and McSaveney, 1983; Reference Henderson and ThompsonHenderson and Thompson 1999). Several glaciers, including Franz Josef, terminate only a few hundred metres above sea level. Ablation at this altitude can reach annual rates of 20 m w.e. (Reference Anderson, Lawson, Owens and GoodsellAnderson and others, 2006), a factor that contributes to a high mass-balance sensitivity.
Until recently, it was not possible to test the idea that New Zealand glaciers are especially sensitive to climatic variations, because few data were available. This paper employs a detailed set of climate and mass-balance data from 2004–08 to construct and test a mass-balance model for Brewster Glacier. The model uses a detailed climatology of the glacier to provide insight into the partitioning of energy-balance components (radiation, turbulent fluxes and precipitation heat) and their influence on the seasonal evolution of the glacier mass-balance profile. The model is used to address the following specific questions regarding the response of Brewster Glacier to climate change, with the overall aim of testing the idea that New Zealand glaciers are highly sensitive:
What are the spatial and temporal patterns of mass balance at Brewster Glacier?
What proportion of annual ablation is caused by different components of the energy balance?
What proportion of annual run-off is caused by seasonal snowmelt, glacier melt and rainfall?
How will energy-balance and run-off characteristics differ under various climatic conditions?
Brewster Glacier
Brewster Glacier covers an area of 2.5 km2 on the western side of the Southern Alps of New Zealand, at 44.05° S, 169.26° E (Fig. 1). The glacier ranges in elevation from 1660 to 2400 m a.s.l. and can be divided into two parts, a large (1.9 km2, 80% of total area), gently sloping (mean 11° gradient) trunk below 2000 m a.s.l. which has a southerly aspect, and a small (0.6 km2), steep (mean 31°) portion above 2000 m a.s.l. which has a southwesterly aspect. The glacier occupies 70% of the catchment that drains to a proglacial stream. While not as well known as Franz Josef Glacier 90 km to the north (Fig. 1), Brewster Glacier has attracted considerable interest. The glacier was mapped roughly by Sir Julius von Haast in 1863 (Reference BurrowsBurrows, 2005), and a photograph taken in 1955 indicates it was then ∼500 m longer than at present. The snowline has been monitored using oblique aerial photography since 1977 (Reference Chinn, Heydenrych and SalingerChinn and others, 2005). Used as a mass-balance proxy, snowline measurements suggest the glacier has experienced positive balance in 16 of the 24 years for which measurements could be obtained (Reference Chinn, Willsman and SalingerChinn and others, 2006). Despite the dominance of positive balance years, the glacier has shown a general retreat of ∼150 m between when the glacier was mapped in 1986 (NZMS 260, G38) and its 2006 position (Fig. 1). This apparently contradictory behaviour, retreating while experiencing positive mass balances, is a result of the long response time of the glacier (estimated at ∼50 years using the method of Reference Jóhannesson, Raymond, Waddington and OerlemansJóhannesson and others, 1989), which means that the glacier is still retreating in response to 20th-century warming, which has been most pronounced since the 1940s (Reference SalingerSalinger and others, 1995).
Fig. 1. Brewster Glacier, and its location in the South Island of New Zealand. Stream stage is measured at S. Climatic variables are measured at the glacier terminus (BG1), at a site at 1760 m a.s.l. (BG2) for rain-gauge measurement, and at a site at 1885 m a.s.l. (BG3) for albedo and snow surface height. Long-term climate data are measured at Haast (H). Ivory (IG) and Franz Josef Glaciers (FJG) are referred to in the text.
Methodology
Climate data
In June 2004 an automatic weather station (designated ‘BG1’) was installed at 1650 m a.s.l., ∼200 m in front of the glacier (Fig. 1). For measuring wind speed an A200m vector anemometer is used, and for wind direction a W200p vector wind vane. Both instruments are mounted 3.5 m above the ground. Air temperature and relative humidity are measured with a SKH 2031 sky temperature and humidity probe, incoming solar radiation is measured with an LI-COR– PY200SA pyranometer and net radiation is recorded with a REBS Q*6.1 net radiometer from Radiation and Energy Balance Systems Inc, all mounted 2 m above the ground. During the winter of 2004, heavy snowfall buried the lower instruments, so in November 2005 the lower sensors were raised to 3 m above the ground. An estimated mean height of these sensors above the ground, taking into account winter snow accumulation, is 2 m. Measurements are taken every 30 s and the climate data are averaged every hour, except for wind speed where an hourly sample measurement is taken.
Precipitation was measured at a storage rain gauge (BG2) located on the western margin of Brewster Glacier at 1760 m a.s.l. (Fig. 1). An Odyssey Data Systems capacitive sensor and data logger records the fluid level in the storage gauge every 15 min. The capacitive sensor is temperature-sensitive, which can lead to false measurements (both positive and negative) during temperature fluctuations. Manual inspection of the rainfall record, in conjunction with the temperature record from the nearby climate station, BG1, enables removal of the majority of the temperature-related errors. Unfortunately, the effect of temperature fluctuations during precipitation events could not be corrected for. The error related to this effect was considered negligible compared to that related to undercatch.
The rain-gauge location was selected to enable precipitation measurement at an elevation as close as possible to the estimated equilibrium-line altitude (ELA) of the glacier while being on stable ground and without the risk of being buried by heavy snowfall or destroyed by avalanche. The exposed location of the rain gauge, the height of its orifice, the lack of wind shield and the frequent occurrence of solid precipitation all lead to the likelihood of significant under-catch. A conservative estimate of the undercatch was obtained through a parameterization, following Reference YangYang and others (1998), based on temperature and wind speed during precipitation events. The temperature at the gauge was estimated by applying a lapse rate to the measurements obtained at BG1. Wind speed was taken to be the same as that obtained at BG1. The estimated undercatch ranges from 59% (high winds and mixed snow and rain) to 9% (light winds, rain), with a mean of 20%.
Climate station BG3, at 1885 m a.s.l. on the glacier surface (Fig. 1), recorded incoming and outgoing shortwave radiation hourly using LI-COR Li-200s and Li-200x sensors facing upwards and downwards 3 m above the snow surface. The data were used to calculate the seasonal evolution of glacier surface albedo. Silicon pyranometers such as these are not recommended for use in measuring reflected radiation, but the measurements are not used as model input, rather for comparison with an established albedo model.
The input data required for the mass-balance and run-off models are temperature, precipitation, wind speed, relative humidity and incoming solar radiation. Annual mean values of these variables are summarized in Table 1. These records do not cover the full 4 year period from 1 April 2004 to 31 March 2008, and the gaps are filled using data from longterm lowland climate stations.
Table 1. Annual climate values measured at climate station BG1 near the terminus of Brewster Glacier. Periods of missing climate data are filled as explained in the text. All quantities are annual averages, except for precipitation which is an annual total
The temperature record was extended using hourly data collected at Haast, on the Tasman Sea coast, 40 km to the west (Fig. 1), using a mean lapse rate of −0.0053°C m−1 calculated between the Haast climate station and BG1. This value is close to the moist adiabat and similar to that found at Franz Josef Glacier, 90 km to the north (Reference Anderson, Lawson, Owens and GoodsellAnderson and others, 2006).
Gaps in the precipitation record for Brewster Glacier were reconstructed using an interpolated precipitation dataset. These data define a two-dimensional surface of daily precipitation interpolated from long-term (and mostly lowland) stations using a trivariate spline, scaled using a subjectively based mean annual precipitation surface to improve estimates at high elevation (Reference Tait, Henderson, Turner and ZhengTait and others, 2006). Comparison of the time series of rainfall for the gridpoint that covers Brewster Glacier with short-term rain-gauge measurements at BG1 indicates that the pattern of daily rainfall is consistent between measured and interpolated values, but that the amount of precipitation is underestimated. As the interpolated surface is calculated daily, rainfall is assumed to fall uniformly through the day to match the hourly time-step of the mass-balance model. The interpolated precipitation value was scaled to the precipitation measured at BG2 from February 2005 to January 2006, using a factor of 1.3.
Analysis of the wind-speed record from BG1 showed that the sensor deteriorated during 2006, recording significantly lower wind speeds in the last 2 years of the study period. Comparison of the BG1 wind-speed record with a record interpolated from nearby long-term climate stations (personal communication from A. Tait, 2008) indicated that during the period when the sensor was operating correctly the mean wind speed at BG1 was 77% of the interpolated value. To ensure consistency between the calibration (2004–06) and evaluation (2006–08) periods the interpolated wind-speed record multiplied by 77% is used for the entire period.
Gaps in the relative-humidity record are filled using uncorrected values measured at Haast. This is unlikely to affect the results significantly, because data gaps occur during winter months when ablation is minimal. As incoming solar radiation was not measured at either BG1 or Haast between April and June 2004, this gap was filled using BG1 data from 1 year later. This approach is preferred over calculations of radiative terms directly as it incorporates the effect of cloud cover.
Energy balance is calculated on a 10 m square grid on Brewster Glacier based on a digital elevation model (DEM) constructed from global positioning system (GPS) data collected in 1997 (Reference Willis, Lawson, Owens, Jacobel and AutridgeWillis and others, 2009) on the glacier surface. The surrounding bedrock topography is based on the 1986 contour map (Fig. 1).
At each time-step, a grid of climate variables over the model domain is created for input to the model. Incoming solar radiation is modified according to shading at each time-step. Precipitation, wind speed and relative humidity are extrapolated uniformly over the grid. A temperature grid is created using a constant lapse rate of −0.0053°C m−1, measured between Haast and Brewster Glacier.
Mass-balance data
To measure ablation, stakes were initially installed on the glacier in February 2004, and extended to a network of 33 stakes (Fig. 1) in December 2004. This network has been maintained to the present. The stakes were set up in a longitudinal profile as recommended by Reference Østrem and BrugmanØstrem and Brugman (1991) and Reference Kaser, Fountain and JanssonKaser and others (2003). The stakes were placed ∼100 m apart with a line of 19 stakes going up the centre of the glacier from 1650 to 1950 m a.s.l. At four locations, transverse profiles of two to four stakes were added along a contour. During summer, the stake network is complemented by a sonic ranger 50 (SR50) that measures the snow height change automatically at BG3 (Fig. 1), next to stake 12.
Accumulation was measured by digging three to four snow pits along the centre line each spring from 2004 to 2007, and one snow-core measurement at 2312 m a.s.l. in 2006. Snow density was measured at the side of the snow-pit walls to calculate accumulation water equivalent. Accumulation was also calculated by measuring stakes from the previous year which reappeared from the snowpack during the ablation season.
The annual mass balance was calculated from ablation-stake, snow-pit and snow-probing measurements. The mass balance was interpolated manually with the help of typical melting patterns revealed by annual oblique end-of-summer snowline photographs since 1977 (e.g. Reference Chinn, Heydenrych and SalingerChinn, 2005) and DEM-derived topographical characteristics such as contour lines, slope, curvature and aspect. The winter balance for the mass-balance years 2005–07 has been interpolated with two linear gradients from snow-pit measurements, applied to the DEM.
Streamflow data
Stream discharge is calculated by measuring stream stage, and relating this to discharge with a rating curve, where discharge is measured independently and related to measured stage. Stage data were collected ∼50 m from the glacier snout (Fig. 1), where the proglacial stream runs through a bedrock channel. An Odyssey capacitive water-level probe mounted in a 2.4 m galvanized-steel pipe was used, recording at 40 min intervals.
A hand-held current meter (type OSS-PC1, serial No. 92-22) was used for flow measurement (Reference Mosley, McKerchar and MaidmentMosley and McKercher, 1993; WMO, 1994). Discharge was measured at ∼1200 and 1800 h, local time, daily over the period 15 February 2006 to 6 March 2006. A discharge vs stage rating curve was created according to the method of Reference Mosley, McKerchar and MaidmentMosley and McKerchar (1993) for the period of direct discharge measurement, allowing interpolation of discharge values over the entire period of stage measurement. There is a very strong relationship (r 2 = 0.95) between discharge and stage. The rating curve is a polynomial, D = 1.335 × 10−5 S 2 − 0.007611S + 1.138, where S is the stage, and is considered to be robust given the stability of the channel at the measurement site (Reference WardWard, 1967; Reference Mosley, McKerchar and MaidmentMosley and McKerchar, 1993; WMO, 1994) and the high proportion of total discharge captured at the measurement site. As a small amount of discharge leaves the main channel before the measurement site, discharge was measured in this subsidiary channel and is included in the rating curve.
Energy-balance model
Melt is calculated using an energy-balance model:
where Q m is the energy available for melt, I the incoming shortwave radiation, α the surface albedo, L out the outgoing longwave radiation, L in the incoming longwave radiation, Q H the sensible heat flux, Q E the latent heat flux, Q R the heat flux supplied by rainfall and Q G the heat flux conducted from the ice. The various components of the energy balance are calculated using climatic and topographic data.
Shortwave radiation
Clear-sky incoming solar radiation, I, is comprised of a direct component, I dir, and a diffuse component, I dif. The direct component is calculated for every unshaded point in the catchment using sinusoidal eccentricity, with the solar constant, S, changing with the Earth–Sun distance (Reference OerlemansOerlemans, 1992):
where ψ a is the atmospheric clear-sky transmissivity, N is the day number and n the cloudiness. The angle of incidence between the slope normal and the solar beam, θ, is calculated as (Reference Garnier and OhmuraGarnier and Ohmura, 1968):
where β is the slope angle, Z the solar zenith angle, φ sun the solar azimuth angle and φ slope the slope azimuth angle (aspect).
The diffuse component of the incoming shortwave radiation, I dif, can be estimated by (Reference OerlemansOerlemans, 1992):
The proportions of direct and indirect solar radiation are then adjusted according to
where the transmissivity of the atmosphere, t a, and of clouds, t c, are given by:
and
where z is the surface elevation. Cloudiness, n, is calculated by an automated trial-and-error approach where the cloudiness is increased from zero until the calculated incoming solar radiation, I, is equal to the measured incoming solar radiation, I m. To avoid having to undertake this calculation at every time-step, a fitted cubic relationship between cloudiness and the fraction of clear-sky to measured radiation is used, a method similar to that described by Reference HockHock (2005), such that:
The constants in Equation (8) are derived from 1 year of hourly measured incoming radiation, I m, at Brewster Glacier with I calculated as above.
The magnitude of the outgoing shortwave radiation depends on the surface albedo, which varies in space and time. Here, we use the method of Reference Oerlemans and KnapOerlemans and Knap (1998) where albedo, α, is calculated as:
where α snow is the albedo of snow, α ice the albedo of ice (0.34), d the snow depth (mm w.e.) and d * the characteristic snow depth scale (11 mm w.e.). The albedo of snow is calculated from:
where α firn is the albedo of firn (0.53), α frsnow is the albedo of fresh snow (0.9), s is the time since the previous snowfall event (days) and t* is a timescale (21.9 days). The values of the constants, taken from Reference Oerlemans and KnapOerlemans and Knap (1998), achieved a good match between modelled and measured albedo at BG3 (see Fig. 1 for location) which remained snow-covered for the period of this study (Fig. 2).
Fig. 2. Measured and modelled albedo at stake 12 at 1885 m a.s.l. The albedo is modelled using the method of Reference Oerlemans and KnapOerlemans and Knap (1998) without modification.
Longwave radiation
Longwave radiation can be calculated for a black-body object using the Stefan–Boltzmann law:
where L is the longwave radiation, ∊ the emissivity of the object, σ the Stefan–Boltzmann constant and T the absolute temperature of the object. We can calculate the outgoing longwave radiation, L out, of the glacier surface, assuming that the surface is at melting point, and the emissivity of snow, ∊ s, is 1. In this case, L out is constant at 317 W m−2.
It is more difficult to calculate incoming longwave radiation, as there are components from the atmosphere and from the surrounding terrain. In this case, the longwave incoming radiation is partitioned according to the viewfield, v, for each gridcell, calculated using the algorithm of Reference CorripioCorripio (2003). Longwave incoming radiation is then calculated as the sum of the atmospheric and terrain components:
where ∊ eff is the effective atmospheric emissivity given by Reference Konzelmann, van de Wal, Greuell, Bintanja, Henneken and Abe-OuchiKonzelmann and others (1994):
and where ∊ c is the clear-sky emissivity given by ∊ c = 0.23 + 0.484(e a/T a)1/8. Here e a is the vapour pressure, T a the air temperature and ∊ c the emissivity for overcast conditions, set at 0.924. The exponent, p = 4, was used by Reference Konzelmann, van de Wal, Greuell, Bintanja, Henneken and Abe-OuchiKonzelmann and others (1994). Measurements of longwave incoming radiation, as the residual from a net all-wave radiometer and shortwave pyranometer for the short period of reliable measurement at Brewster Glacier, indicate that p = 4 tends to underestimate longwave incoming radiation, so we set p = 1. The terrain emissivity, ∊ t, is set to 0.4 (Reference Plummer and PhilipsPlummer and Philips, 2003). The temperature of the terrain, T t, is taken as the atmospheric temperature (corrected for elevation) for terrain that is not snow- or ice-covered, or 273 K for terrain which is snow- or ice-covered.
Turbulent heat fluxes
The turbulent fluxes of sensible and latent heat are calculated using the bulk aerodynamic approach which exploits the well-defined conditions of a melting surface (Reference OerlemansOerlemans, 1992):
where ρ is the density, c p the specific heat capacity of ambient air, k H and k E the exchange coefficients, and U the wind speed at height z = 2 m. The temperature of the surface, T s, is assumed to be 0°C, L v is the latent heat of vaporization for water, q the vapour pressure of the ambient air, q s the vapour pressure of air at the glacier surface and p the air pressure. The exchange coefficient for sensible heat (and similarly for latent heat) is related to the roughness length for wind, z 0, by:
where k 0 = 0.4 is von Kármáńs constant, z 0H the roughness length for sensible heat and R b the Richardson stability criterion:
where g is acceleration due to gravity (9.8 m s−2).
The stability correction (the final term in Equation (16)) is only applied when R b > 0 and when U > 1, the latter to prevent the bulk Richardson number from becoming unreasonably high when wind speeds are very low (Reference OkeOke, 1987). This is a reasonable approach given that at these times the turbulent heat fluxes will also be low. It is assumed that the roughness lengths for wind, sensible and latent heat fluxes are equal. We treat this value, z 0, as an ‘effective roughness length’ and tune it to provide appropriate values for the exchange coefficients (Reference BraithwaiteBraithwaite, 1995).
Rainfall and ground heat fluxes
The precipitation heat flux is given by:
where c w is the specific heat of water and P the rainfall rate. Rainfall is assumed to be at air temperature, T a.
The ground or ice heat flux, Q G, can be important, especially in polar environments, but requires calculation of subsurface energy flows. We assume that the surface stays at 0°C. While not strictly correct, this is expected to have little influence on the mass-balance calculation on temperate glaciers (Reference OerlemansOerlemans, 1992). The ground heat flux is set at Q G = 1 W m−2 (Reference Neale and FitzharrisNeale and Fitzharris, 1997).
Accumulation model
Accumulation is modelled using a simple temperature threshold to discriminate rain from snow. The threshold, T s, is set at 1°C, a value used in other snow-modelling studies in New Zealand (Reference Moore and OwensMoore and Owens, 1984; Reference BarringerBarringer, 1989; Reference Anderson, Lawson, Owens and GoodsellAnderson and others, 2006). The sensitivity of model results to this choice is tested below. Application of this model to areas off the glacier resulted in unrealistically deep ‘seasonal snow’ (as indicated by photographs of late-winter snow cover). Without measurements of seasonal snow accumulation a simple scheme, where snow accumulation is halved in areas off the glacier, was adopted. This adjustment only affects the run-off calculation and will be revisited in that context.
Run-off model
Discharge is calculated using a linear reservoir model (Reference Baker, Escher-Vetter, Moser, Oerter and ReinwarthBaker and others, 1982; Reference Hock and NoetzliHock and Noetzli, 1997). This is a ‘lumped’ model, where the water inputs are not spatially distributed. Despite this limitation, the model often performs as well as complex models (Reference HockHock, 2005).
In this model, the discharge, D(t), at time t is proportional to the reservoir’s volume, V(t):
where k s is the storage constant. The rate of change of volume is given by:
where R(t) is the rate of water inflow into the reservoir. Combining these equations and solving for D, we have:
The glacier is split into three reservoirs, each with a different storage constant, k s, one for snow, one for firn and one for ice. Values for the storage constants are taken from Reference Hock and NoetzliHock and Noetzli (1997) with k snow = 350, k firn = 30 and k ice = 16.
Model tuning and evaluation
The model is tuned by comparing modelled mass balance against measured mass balance. The aim is to tune as few parameters as possible. As we have sufficient measurements to constrain the albedo calculation and the temperature lapse rate, these are untuned. The snow/rain threshold is unmeasured but we use accepted values so have not tuned this value. The sensitivity of the model to the tuned and untuned parameters is investigated.
Energy-balance model tuning
Measured and modelled mass balances are compared for 724 individual mass-balance measurements over the first 2 years (April 2004–March 2006) of the 4 year study period (Fig. 3). Local mass balance is best known at stake 12 (1885 m a.s.l.), where the sonic ranger was installed in late summer in 2005 and 2006. This stake is near the steady-state ELA (as defined by Reference Chinn, Heydenrych and SalingerChinn and others, 2005) at 1935 m a.s.l. (personal communication from T.J. Chinn, 2008). A time series for this stake (Fig. 4) indicates the mass balance is well simulated through this period. The match between modelled and measured mass balance has been achieved by adjusting two parameters.
Fig. 3. Measured and modelled mass balance compared for all mass-balance measurements. Best-fit lines have been calculated separately for the calibration period, April 2004–March 2006 (r 2= 0.76, root-mean-square error = 0.03), and the evaluation period, April 2006–March 2008 (r 2= 0.72, root-mean-square error = 0.40).
Fig. 4. Mass-balance model output and measured mass balance at stake 12 at 1885 m a.s.l. At the start of the record and at measurement gaps, the measured mass-balance record is aligned to the modelled mass-balance record and marked ‘+’. The sonic ranger data are converted to mass-balance changes using a constant density of 600 kg m−3.
Precipitation, measured at 1760 m on the glacier margin, has been adjusted by multiplying it by a factor of 1.2, reflecting an underestimate of accumulation. This factor provides an adequate fit for both short-term accumulation events as recorded by the sonic ranger at stake 12 (Fig. 4) and long-term accumulation over two winters. The factor may represent undercatch at the rain gauge, despite the Reference YangYang and others (1998) correction undertaken (described above), or influences on accumulation such as wind redeposition of snow or avalanching.
A single roughness length for wind, z 0, could not accurately simulate turbulent heat flux on both the snow-and ice-covered portions of the glacier. Tuning experiments resulted in values of z 0 = 0.008 m for snow surfaces, and z 0 = 0.012 m for ice surfaces. These values are towards the upper end of measurements on snow and ice surfaces (Reference Brock, Willis and SharpBrock and others, 2006).
Run-off model evaluation
A time series of discharge has been calculated using the stage data and rating curve described above (Fig. 5). While the model simulates the temporal patterns of discharge, and at times the magnitude of rated discharge, closely, in general the rated discharge is lower than the modelled discharge. The most likely reason for this result is that parts of the catchment where accumulation and ablation have not been measured are behaving differently in response to the low-angle glacier surface. In particular, the off-glacier seasonal snow accumulation has not been measured and its melt can have a significant impact on summer discharge. The off-glacier seasonal snow accumulation was assumed to be equal to half that observed on the glacier. To test the sensitivity of this assumption, model runs were made with the off-glacier seasonal snow set to zero and also set to the same value as the on-glacier accumulation. With no off-glacier snow accumulation, the overestimation of the modelled discharge dropped from 55% to 34%. With the accumulation levels the same on and off the glacier, the modelled discharge increased to 69% more than the observed discharge. These results indicate the off-glacier accumulation estimation partially explains the observation-to-model mismatch, but as a zero off-glacier accumulation is unrealistic, there are clearly other effects unaccounted for. The remaining difference between modelled and measured discharge shows some variation through the melt season; the magnitude is approximately correct up until late January but is overestimated by an increasing amount through the summer. There are insufficient data to further constrain the reasons for this difference.
Fig. 5. Modelled, measured and rated discharge for mid- and late summer 2006. The range of measured discharge is small, and extrapolation beyond this range is necessary for higher discharges. In the period 22 January to 13 March 2006 the mean rated discharge was 0.81 m3 s−1, 64% of the mean modelled discharge of 1.26 m3 s−1.
Model sensitivity
The parameters used in tuning were systematically varied to assess the robustness of the model. The sensitivity to the temperature lapse rate and albedo calculation were also assessed, although these have not been tuned (Table 2).
Table 2. Deviation of mass balance from standard run for different model parameter values
The results show that the model is sensitive to the rain/snow temperature threshold, T s, surface albedo (snow and ice) specification and the temperature lapse rate, dT/dh. The sensitivity of the model to these parameters is a reflection of the high precipitation rates, where a small change in how precipitation is partitioned between rain and snow can lead to a large change in the magnitude of accumulation experienced by the glacier.
Results
Patterns of mass balance
Local annual mass balance decreases from approximately +3 m w.e.a−1 near the highest point on the glacier to approximately −1 to −4 m w.e.a−1 at the terminus (Fig. 6). There is little ablation at the highest elevations and most precipitation falls as snow. On the gently sloping glacier tongue, there are variations in mass balance which are independent of elevation. For example, mass balance is less negative on the eastern margin than on the western margin at the same elevation. The pattern also shows up in oblique aerial photographs of Brewster Glacier taken at the end of the summer (Fig. 7), and is the result of increased shading and decreased ablation on the eastern margin.
Fig. 6. Contour map of mass balance (m w.e. a−1)modelled over the glacier surface for each year, finishing 31 March. The spatial variations of mass balance are very similar for the four years, but the magnitude of both accumulation and ablation (and hence the mass-balance gradient) change significantly from year to year.
Fig. 7. Brewster Glacier at the end of summer 1998. The eastern margin is to the right of the photograph. Snow- and firn lines from three previous years, but influenced by melting in 1998, are visible. The shape of the lines indicates that the mass balance experienced by the glacier varies in a complex way. (Source: T.J. Chinn.)
Annual energy-balance components
Summing the energy-balance components at monthly resolution allows the sources of energy and their seasonal evolution to be identified (Fig. 8). In the winter months (June, July and August) small amounts of melt are caused almost exclusively by radiation (Fig. 8), as the air temperature is too low to allow much, if any, energy transfer to the surface by turbulent fluxes. Turbulent transfers are the primary source of energy for melt during the summer. Precipitation also contributes to melt in the summer, as the portion of precipitation falling as rain increases, and the temperature of the rain (assumed to be the same as air temperature) increases.
Fig. 8. Variation in source of energy for melt throughout the 4 year study period. Asymmetry in the annual pattern results from glacier surface feedbacks. The abrupt transition from dominantly ice to snow, which occurs in late summer or autumn, causes an increase in albedo and a decrease in surface roughness over large parts of the glacier simultaneously, resulting in reduced energy fluxes and hence lower ablation.
On an annual scale, slightly more than half (mean over 4 years of 52%) of the energy available for melt comes from turbulent heat fluxes (Table 3) and radiation provides 45% of melt energy. Even in this high-precipitation environment, precipitation contributes only a few per cent of the total heat flux.
Table 3. Annual energy-balance components for each of the four years. The proportion of the total for each component is shown in parentheses
Annual run-off components
Run-off is separated into three components: glacier melt (defined as the melt of snow and ice inside the glacier boundary), seasonal snowmelt (defined as the melt of snow outside the glacier boundary) and rainfall (over the entire catchment). Each component responds differently to the climate, and involves storage of water on different timescales.
Over the 4 year period the model shows a strong seasonality of run-off, with June to August in each year having a discharge close to zero (Fig. 9). Rapid increases in discharge during spring are due mostly to rainfall run-off, as the air temperature rises and a greater proportion of precipitation falls as rain rather than snow. Snow- and ice melt during the spring increases more slowly, reaching a peak between December and February. Towards the end of the summer, the contribution of seasonal snow declines as surrounding bedrock topography is exposed. Glacier melt continues to be an important source of meltwater into March, until the first significant snowfalls in April and May effectively terminate the summer melt cycle.
Fig. 9. Components of melt and run-off calculated by the energy-balance and discharge models. For clarity a 1 month running mean is used to smooth the discharge curve.
On an annual scale, rainfall contributes more to run-off than glacier melt, despite the fact that the catchment is 70% glacierized. The contribution of seasonal snow to run-off is a distant third, reflecting both the smaller proportion of the catchment that is not glacierized and the seasonal snow lost from the catchment by the end of summer.
Mass-balance sensitivity to climate change
The Brewster Glacier total mass balance, that is the mass balance summed over the glacier surface for a balance year, is extremely sensitive to temperature change (Fig. 10), especially on the 80% of the glacier below 2000 m (Fig. 1). For Brewster Glacier, the climatic conditions in 2004–08 indicate a mean temperature sensitivity for the 4 years of −2.0 m w.e.°C−1. This is calculated using the mean of the mass-balance change for +1 and −1°C (Reference Braithwaite, Zhang and RaperBraithwaite and others, 2002). The extremely high sensitivity is a result of the limited altitudinal range of the majority of the glacier (500 m; Fig. 1) and glacier surface feedbacks. As temperature increases, higher turbulent heat fluxes lead to higher ablation rates, which leave more of the glacier as exposed ice. Exposed ice absorbs more radiation because of its lower albedo, and the turbulent fluxes increase because of higher surface roughness.
Fig. 10. The sensitivity of annual mass balance to (a) temperature and (b) precipitation changes for the 4 years 2004–08. The dramatic response to temperature change, evident below 2000 m a.s.l., indicates that the glacier is unlikely to survive a modest warming.
For precipitation, the mean mass-balance sensitivity for the climate conditions in 2004–08 is 0.4 m w.e. for a 10% change. Hence a 50% change in precipitation results in the same change in total mass balance as a 1°C change in temperature.
Run-off sensitivity to climate change
Run-off from the glacier is also extremely sensitive to temperature change. A 1°C increase in temperature leads to a 60% increase in run-off; conversely a 1°C decrease in temperature leads to a 20% decrease in run-off. The high sensitivity of run-off to higher temperatures is also the result of feedbacks associated with albedo and turbulent fluxes. The lower albedo and higher surface roughness of an ice surface means that ablation rates increase as higher temperatures change the proportion of snow-covered to ice-covered areas on the glacier surface.
For precipitation, a 20% increase in precipitation leads to only an 11% increase in run-off, while a 20% reduction in precipitation leads to only a 6% decrease in run-off. In this case, there is a strong negative feedback, whereby an increase in precipitation increases the amount of snow and ice stored in the catchment (positive mass balance). A decrease in precipitation leads to a decrease in snow and ice storage, replacing run-off that would have come directly from rainfall.
Discussion
Spatial and temporal patterns of mass balance
As on most alpine glaciers, the dominant variable influencing local mass balance is elevation. This variation is predominantly a result of the reduction in temperature with elevation, which influences the turbulent transfers, incoming longwave radiation and partitioning of precipitation between snow and rain. The use of a distributed model allows assessment of mass-balance variations that are independent of elevation. In this model the radiative energy fluxes are independent of elevation and are influenced by other terrain factors.
Shading from surrounding topography and the glacier slope and aspect influence the magnitude of direct incoming shortwave radiation, and the sky-view factor influences the amount of diffuse radiation received. The effect of terrain on longwave radiation is less important, due only to the sky-view factor, which changes the relative amount of radiation received from the atmosphere and surrounding terrain. On other glaciers shading and the combined effect of aspect and slope share equally in the effect on energy fluxes, and are several times larger than the effect of sky view (Reference Arnold, Rees, Hodson and KohlerArnold and others, 2006).
Due to the relatively large (∼52%) proportion of energy supplied by turbulent heat fluxes, modelled patterns of mass balance generally strongly follow elevation (Fig. 6), with the exception of an east–west gradient in mass balance on the glacier toe.
The pattern of mass balance implied by oblique aerial photographs of the glacier at the end of summer (Fig. 7) appears to be more complex than the model predicts. Some of this variability may result from simplified modelled ablation patterns resulting from the smoothed nature of the DEM (cf. Reference Arnold, Rees, Hodson and KohlerArnold and others, 2006) and more complex wind patterns than the uniform pattern assumed. However, there is also greater variability in accumulation than can be modelled by elevation alone.
The model does not take into account snow redistribution from wind or from avalanching, which occurs on the steep slopes around Mount Brewster. The patterns of mass balance on the steep slopes of Mount Brewster have not been measured, so the magnitude or effect of snow redistribution cannot be quantified. To fully understand the patterns of mass balance at Brewster Glacier these processes need to be examined in more detail.
The modelled annual mass balance of Brewster Glacier in the four balance years varied from 0.6 to −1.9 m w.e. a−1 (Table 4), with the overall balance for the 4 years being −0.7 m w.e. a−1. Modelled winter balances (April–October) were consistently within the range 2.3– 2.8 m w.e. for the 4 years, while summer balances (November–March) show a much greater variation, in the range −1.4 to −4.5 m w.e. (Table 4). For these 4 years, at least, we can infer that the interannual variability in total mass balance is driven primarily by the summer balance. In turn, summer balance is driven primarily by temperature and we conclude that, at least for the study period, mass-balance changes were driven by the temperature variations evident in Table 1.
Table 4. Annual mass balance (m w.e. a−1) calculated for Brewster Glacier using the energy-balance model for the four study years. For comparison, total mass balance interpolated directly from stake measurements is also shown (Reference Haeberli, Hoelzle and ZempWGMS, 2007)
There is some variation between total measured mass balance on the glacier, calculated using manual interpolation of the mass-balance measurements, and total modelled mass balance (Table 4). The difference is ∼0.5 m w.e.a−1 in all years, but modelled values are lower than measured in 2004–06 and higher than measured in 2006–08. The differences are a combination of the model inaccuracies (Fig. 3) and observation interpolation inaccuracies, primarily in the winter balance calculations.
Annual energy-balance components and ablation
In general, the radiative energy fluxes are less important as a source of energy for melt in maritime climates than in less cloudy continental climates. Reference Willis, Arnold and BrockWillis and others (2002) review a number of energy-balance studies and categorize them into maritime and continental climates. In general, the radiation flux dominates. In continental areas, radiation comprises 77% of the fluxes, while in maritime areas (including New Zealand) energy for melt is evenly split between the radiative and turbulent fluxes.
Previous energy-balance studies on glaciers and snow surfaces in New Zealand come from point measurements, and indicate that turbulent fluxes generally dominate by a large margin, more so than for other maritime areas. Reference Moore and OwensMoore and Owens (1984) and Reference Moore and ProwseMoore and Prowse (1984) summarize winter and spring energy-balance calculations in nonglacierized catchments to the east of the main divide, and conclude that radiation provides only a small proportion (∼20%) of total energy for melt. Short-term energy-balance measurements during fine periods on the lower Franz Josef Glacier have also found that turbulent heat fluxes dominate in both summer and winter (Reference Ishikawa, Owens and SturmanIshikawa and others, 1992).
In contrast, other high-elevation energy-balance studies in New Zealand (Reference Kelliher, Owens, Sturman, Myers, Hunt and McSevenyKelliher and others, 1996; Reference Neale and FitzharrisNeale and Fitzharris, 1997) suggest that, even in the maritime climate of the Southern Alps, lower temperatures and reduced atmospheric absorption at high elevations result in net radiation being the most important energy source. However, Reference CutlerCutler (2002) found turbulent fluxes contributed 67% of melt energy at a high-elevation site, indicating that the relationship is not simple.
The results of this study need to be treated with caution, as the energy fluxes are not measured but are calculated to fit measured mass balance. However, the results indicate that there are both spatial and temporal variations in the proportion of energy available for melt from different sources. Radiation dominates the energy balance in winter, while turbulent fluxes dominate both in summer, when temperatures are higher, and on an annual scale.
Annual run-off components
Run-off directly from rainfall comprises ∼50% of the annual discharge from Brewster Glacier, at least in the 4 years of this study. This proportion is rather lower than that measured at Ivory Glacier, 120 km to the north (Fig. 1). During 2 years of measurements at Ivory Glacier, the run-off from rainfall was 79% of a total of ∼10–11 m of annual run-off in 1971–75 (Reference AndertonAnderton, 1976). Ivory Glacier is in a significantly wetter and somewhat warmer (lower-elevation) environment than Brewster Glacier.
The proportion of direct run-off from rainfall is much higher for these New Zealand glaciers than for glaciers of the Europeans Alps. For example, Aletschgletscher, Switzerland, which has the same proportion of the catchment glacierized as Brewster Glacier (∼70%), only 15% of the annual run-off is directly from rain (Reference Schaper, Martinec and SeidelSchaper and others, 1999).
Significant quantities of water are stored as snow and ice during the winter and released during the summer at Brewster Glacier. As in other glacierized basins (e.g. Reference Nienow, Sharp and WillisNienow and others, 1998), the seasonality of run-off is a direct result of this storage, as there is little seasonal change in the quantity of precipitation. During the winter months, almost all precipitation is stored as snow, with almost no run-off during June to August. Approximately half the precipitation that falls on the glacier is stored and released in an annual cycle (Table 5). Almost all seasonal snow storage is released by the end of the summer, and ∼31% of the annual rainfall for 2004–08 went into long-term storage as firn and eventually glacier ice.
Table 5. Annual run-off components for each of the four years. The proportion of the total for each component is shown in parentheses
Climate sensitivity
Mass balance
The glaciers in New Zealand experience some of highest precipitation rates and the highest mass-balance sensitivities in the world. The 4 year average climate sensitivity of the mass balance of Brewster Glacier to changes in temperature is −2.0 m w.e.a−1 °C−1. As this value is only calculated over 4 years, it may not be representative of the sensitivity of the glacier under a range of climatic conditions. Nevertheless, we conclude that the sensitivity is extremely high.
Reference Braithwaite, Zhang and RaperBraithwaite and others (2002) report a variation of sensitivity from −0.1 to −1.2 m w.e. a−1°C−1 for 61 glaciers and ice caps and note that higher sensitivity was associated with warm wet maritime environments. The high sensitivity calculated in this study supports that relationship.
Reference OerlemansOerlemans (2001) identified three factors that relate to the increasing sensitivity of mass balance in wet environments: (1) stronger albedo feedback; (2) larger effect on the partitioning of precipitation between snow and rain; and (3) longer ablation seasons because glaciers extend to lower elevations.
There is a strong albedo feedback at Brewster Glacier: changes in temperature and precipitation, through their influence on snow cover and hence albedo, influence the magnitude of outgoing shortwave radiation. For temperature, the total energy flux from net radiation increases by 35% for a 1°C increase (calculated from the ‘No change’ and ‘1°C’ cases in Table 6). A 20% reduction in precipitation results in an increase in radiation energy flux of 14% due to the albedo feedback.
Table 6. Energy-balance components over the period 2004–08 under different imposed climate changes. The proportion of the total for each component is shown in parentheses
It is clear there is a second feedback of a similar magnitude operating at Brewster Glacier, in the turbulent heat fluxes. The turbulent heat fluxes increase at almost the same rate as net radiation heat fluxes for warming: a 1°C warming causes a 33% increase (Table 6). Of course, the magnitude of the turbulent heat fluxes is directly influenced by temperature, but there is also an indirect feedback as ice has a higher roughness than snow and more of the glacier is ice-covered in warmer or drier conditions. This feedback is more important in a maritime climate as there is already a high proportion of total energy available for melt coming from turbulent heat fluxes. For Brewster Glacier (for example) the turbulent heat increase is of the same order as the increase in net shortwave radiation.
Insight into the importance of the partitioning of precipitation between rain and snow may be obtained by comparing the effect of imposing a temperature change of +1°C on the model with changing the snow/rain threshold temperature to T s = 0°C (instead of the model reference value of 1°C). The mean mass-balance deviation over the 4 years of study resulting from changing the threshold temperature is −1.2 m w.e.a−1 (Table 2). This value is 60% of the temperature sensitivity of the model, indicating that the partitioning of precipitation into rain and snow, with associated albedo and turbulent heat feedbacks, is responsible for the bulk of the change in mass balance due to warming.
The final reason for the high sensitivity of glaciers in wet climates put forward by Reference OerlemansOerlemans (2001) is longer ablation seasons because glaciers extend to lower elevations. This is not an important factor at Brewster Glacier, because the glacier terminates at ∼1660 m a.s.l., unlike, for example, Franz Josef Glacier which experiences year-round ablation on its tongue (∼300 m a.s.l.).
However, Brewster Glacier has a geometry that makes it more sensitive to climatic changes than Franz Josef Glacier. The low surface gradient of the glacier means that ∼80% of the glacier area is vertically within 300 m of the terminus (Fig. 1). Relatively minor changes in temperature and precipitation can change this area from a predominantly snow surface in the summer to an ice surface (Fig. 10). This limited elevation range also makes the partitioning of precipitation between rain and snow over the glacier surface sensitive to changing atmospheric temperature.
Our findings provide insight into the major drivers of recent mass-balance changes on New Zealand glaciers. It is known that changes in atmospheric circulation associated with the El Niño–Southern Oscillation influence glacier mass balance (Reference Hooker and FitzharrisHooker and Fitzharris, 1999). This is evident in annual ELA measurements carried out since 1977 (Reference Chinn, Willsman and SalingerChinn and others, 2006). During the study period, El Niño conditions were experienced in the 2004/05 and 2006/07 balance years, with a weak La Niña in 2005/06 and a strong La Niña in 2007/08. El Niño atmospheric circulation in New Zealand generates anomalous southwesterly airflow which can result in higher precipitation, lower temperatures and increased cloudiness. While temperatures were lower, and precipitation higher in the two El Niño years, the relatively small sensitivity of mass balance to precipitation changes identified in this study indicates that higher precipitation is unlikely to be the dominant variable influencing recent mass-balance trends. It is more likely that the reduced temperature dominates.
Run-off
Changes in total precipitation, including rain and snow, have little influence on annual run-off, because increasing precipitation is added to long-term snow and ice storage, and decreasing precipitation leads to water being released from storage, thus dampening the effect.
The run-off sensitivity described herein neglects the influence of glacier dynamics on mass balance. Glacier dynamics have a significant influence on run-off at longer timescales as the geometry readjusts to a new climatic state: an initial increase in run-off following warming is followed by a decrease as the glacier retreats (Reference Ye, Ding, Liu and LiuYe and others, 2003; Reference Jóhannesson, Aðalgeirsdóttir, Björnsson, Pálsson, Sigurðsson and JärvetJóhannesson and others, 2004; Reference Horton, Schaefli, Mezghani, Hingray and MusyHorton and others, 2006). A coupled mass-balance/glacier-flow model would be required to predict run-off changes more than a decade or so into the future.
Conclusion
Previous energy-balance studies in New Zealand have reported that an extremely high proportion (up to 84%) of energy available for melt comes from turbulent heat fluxes, in contrast to more continental areas where net radiation dominates the energy balance. Our study supports the important role of turbulent heat fluxes in a maritime environment. We also suggest that, while the relative importance of turbulent heat fluxes over radiation diminishes with increasing elevation, the patterns are complex and not amenable to a simple maritime/continental division.
On a global scale, this study confirms that the most maritime glaciers should be the most sensitive to a warming world, and will hence make the largest contribution to sea-level rise in the short to medium term. Glaciers in the high-precipitation environment of the western Southern Alps are extremely sensitive to temperature changes. The dominant effects are albedo and surface roughness feedbacks, the partitioning of precipitation into rain and snow, and the relatively high proportion of melt that is caused by turbulent transfers. These feedbacks are predominantly controlled by temperature. The high sensitivity not only includes glaciers such as Franz Josef Glacier, which show dramatic advance and retreat for relatively small climatic changes, but also the smaller, less dynamic glaciers such as Brewster Glacier. These glaciers have static sensitivities to climate warming that are some of the highest reported in the literature. Such a large sensitivity means that small interannual changes in temperature are likely to drive a large change in glacier mass balance in New Zealand, whereas variations in precipitation are less important.
Sensitivity of run-off to temperature changes is also large in this high-precipitation environment. In the short term, a relatively minor warming will lead to significant increases in run-off in glacierized catchments. In the longer term, run-off will decrease as it is limited by the total ice volume in the Southern Alps. This has implications for the management of hydroelectric power-generation systems in the South Island. Present-day inflows into these systems are elevated due to the decrease in ice in the Southern Alps (Reference Hoelzle, Chinn, Stumm, Paul and HaeberliHoelzle and others, 2007) resulting from a 20th-century warming. As there is a finite store of water in these glaciers, glacier run-off will decrease in the long term.
Acknowledgements
B.A. was supported by Victoria University’s University Research Fund (URF 39180). Fieldwork for B.A. and A.M. was supported by a Foundation of Research, Science and Technology (ANZICE, VICX0704) grant and the Comer Science and Education Foundation. Climate station BG1 was established by the University of Otago and we acknowledge the work of N. Cullen in maintaining the station. D.S. acknowledges the Humanities Division, the Active Earth Processes Research Theme and the Geography Department of the University of Otago for financial support, and thanks the World Glacier Monitoring Service and the 3G group of the Geography Department of the University of Zürich for advice and hospitality during several research visits to Zürich.
