Glacier mass balance modelling

The model calculates the surface mass balance of glaciers over a 38 × 40 km domain in the central Southern Alps of New Zealand (Fig. 3 and Supplementary Fig. 1), at 100 m horizontal resolution. Ablation (all processes leading to ice loss at the glacier surface, dominantly melt but also sublimation) is calculated using a spatially distributed energy balance model10. The model is run at a daily time step, and the results are collated for each glacier year, which starts on 1 April and ends on 31 March of the following year (end of the austral summer). The study covers a period of 39 glacier years from March 1972 to March 2011. The model is forced by daily climate data from an ∼5 × 5 km (0.05° lat/lon) climate product derived from instrumental station data and climate re-analysis data.

Glacier surface mass balance depends not only on climatic variations but also on the dynamic adjustment of glacier geometry. For this reason, reference-surface balance calculated over a constant glacier geometry is better suited for climatic interpretation56. Thus, our glacier model considers the surface mass balance (and not the flow dynamics) of glaciers in the central part of the Southern Alps31, (Fig. 3) over a fixed geometry. A 15 m digital elevation model for New Zealand57 is resampled to 100 m to provide the surface elevation throughout the model domain. An ice mask and glacier debris cover are derived from aerial photography58 captured for topographic mapping (NZ map series 260), and resampled on the same 100-m grid as the digital elevation model.

Ablation is calculated using an energy balance model31. Incoming short-wave radiation is derived using a diffuse and direct component59, from calculated radiation at the top of the atmosphere60. A cloudiness value is inferred from measured incoming solar radiation from the New Zealand Virtual Climate Station Network data set28,29 to partition direct and diffuse components. The spatial distribution of short-wave radiation is calculated using topographic shading61, and surface albedo depends on the length of time since the last snowfall62. Outgoing long-wave radiation (L OUT ) is calculated using the Stefan–Boltzmann law, while incoming long-wave radiation is split into an atmospheric and terrain component according to the sky view of each point on the grid. The emissivity and emission temperature of these components are different, with the atmospheric component based on air temperature, with an emissivity that incorporates assumptions about the vertical temperature profile.

Turbulent heat fluxes, Q H and Q C , are calculated using a bulk aerodynamic approach, with different effective roughness length for snow and ice surfaces (Supplementary Table 1). These values were derived by tuning these effective roughness lengths until melt rates achieved a good match with 455 individual glacier mass balance measurements made on Franz Josef Glacier between 1 April 2000 and 31 March 2002. The Richardson stability criterion is used to correct the exchange coefficients for stable conditions. Precipitation heat flux, Q R , is calculated assuming that precipitation is at air temperature. Subsurface heat fluxes are neglected which, while not strictly correct, is expected to have a relatively minor (∼10%) influence on the mass balance calculation on temperate glaciers, which do not contain significant amounts of snow and ice below 0 °C (ref. 59). The ground heat flux, Q G, is set at 0 Wm−2, consistent with the assumption of surface temperature at 0 °C (ref. 63). Surface debris cover is dealt with in a simple way where ablation calculated by the energy balance model is reduced by 90% where surface debris is present and when any snow accumulation on top of the debris has melted, based on observed reductions in ablation under debris on Tasman Glacier of 93% (ref. 64) and 89% (ref. 65).

To ensure that our results are not influenced by these simplifying assumptions regarding subsurface heat fluxes and melt under debris cover, the experiment is repeated with subsurface thermal calculations in snow, ice, debris and snow on debris. The evolution of temperature T in the snow/ice/debris with time t and at depth z is calculated by solving the one-dimensional heat equation (1):

where ρ is the density of the material, c its specific heat capacity and k its conductivity (values summarized in Supplementary Table 1). The heat equation is solved independently at each grid cell using a Crank–Nicolson scheme and a standard tridiagonal solver. To resolve temperature correctly, the time step is reduced from daily (in the standard energy balance calculation) to 4 hourly to capture the diurnal evolution of the thermal state caused by changing energy fluxes at the surface. Air temperature is allowed to vary sinusoidally within the maximum and minimum temperatures of each day, while the short-wave energy components are calculated explicitly at each time step. Surface temperature is then calculated by closing the energy balance at the surface66 using a Newton–Raphson iterative scheme67 where the energy balance components that depend on surface temperature (Q H , Q E , L out and Q R ) are recalculated to estimate the derivative dQ s /dT s of surface energy Q s with surface temperature T s . The boundary conditions for the temperature calculations in each case are the same, with the surface temperature applied at the top and 0 °C applied at the bottom layer. In the case of snow and ice, we are assuming that the glacier is temperate and that at 10 m below the surface there is no annual variation in temperature.

Debris thickness is poorly known in the study area, and so we take the data from Tasman Glacier68 and extract a relationship between debris thickness and the distance from the closest clean-ice surface. This relationship is then applied to Tasman Glacier and all other debris cover in the model domain. The calculation of conductive heat flux within the debris layer generally follows ref. 67. We assume that the ice below the debris is always at melting point. This assumption is justified during the summer67 when surface temperatures are generally above 0 °C. In thin debris cover during the winter the base of the debris will drop below 0 °C, but the amount of ablation under debris cover during these conditions is zero. We also assume that the base of any snow that accumulates on debris has a temperature of 0 °C, which also suppresses any ice ablation at the base of the debris layer. Snow on debris has no direct effect on glacier mass balance (debris cover is always in the ablation zone, so snow on debris never accumulates from year to year) but the duration of snow cover controls melt of the ice underlying the debris. Snow on debris is rather ephemeral and so this treatment is considered adequate. Ultimately, melt under debris is calculated based on the conductive heat flux at the ice/debris interface.

The data used as input to the model are minimum and maximum air temperature, solar radiation, relative humidity and precipitation. As wind speed data are not available for the whole period, wind speed at the 850-hPa level from climate re-analysis30 is scaled to match wind speed from the New Zealand Virtual Climate Station Network for 2000–2010, and applied throughout the model period. All of the climate data fields are downscaled using bilinear interpolation from the 0.05° New Zealand Virtual Climate Station Network grid to a 100 m resolution square grid for the glacier mass balance model. Minimum and maximum temperatures are calculated from lowland stations using a seasonally variable lapse rate, which is different for minimum and maximum temperatures28,69. To interpolate these data to the finer glacier grid, the temperature values are first lapsed to the sea level using the same lapse rates before bilinear interpolation, and then lapsed to the 100-m grid elevations. In this way the dependence of temperature on elevation is captured, besides the regional variations in temperature.

We also carry out an experiment where the mass balance model is forced by temperature data downscaled from climate re-analysis data30 (Fig. 6 and Supplementary Fig. 4), rather than directly from the New Zealand Virtual Climate Station Network. This allows us to assess the utility of using climate reanalysis data to identify the climate drivers of glacier advance and retreat. Downscaling is based on diagnostic regression relationships between surface temperature over New Zealand and variables taken from reanalysis data at 2.5° latitude/longitude resolution using the approach described in ref. 70. Regression relationships are developed for each standard 3-month season for surface temperature, using re-analysis 1,000 hPa height and 850 hPa temperature fields. For development of the regression equations, temperature data are taken as anomalies from daily climatological values. For maximum, minimum and mean temperature, the daily climatology accounts for up to 70% of the total variance, while the re-analysis-based regression estimates account for between 40 and 60% of the residual in western and alpine regions. Downscaled temperatures derived from re-analyses are then extrapolated to our finer-scale grid for glacier mass balance modelling, using the methods described above.

Model evaluation and uncertainties

To test the veracity of the climate date interpolation schemes, we compare the model input data against independently measured data at Mueller Hut and Tasman Glacier. We limit this comparison to temperature because there are very few reliable precipitation measurements in the model domain, and accumulation measurements (tested later) are the most robust way of assessing precipitation input. The Tasman Glacier temperature data71 were collected from 2007 to 2009 at 1,139 m above the sea level, and have a root-mean square difference of 2.32 K and a bias of 0.86 K compared with our model input data. The Mueller Hut temperatures72 were collected from 2010 to 2011 and show root-mean square difference of 2.17 K and a bias of 0.30 K. These statistics compare favourably with other temperature interpolation methods in the Southern Alps28.

Glacier mass balance measurements made at Franz Josef and Tasman glaciers are used to tune and evaluate the glacier mass balance model. The model is also evaluated against annual glacier snowlines in the Southern Alps (a proxy for the glacier equilibrium line altitude) estimated from oblique aerial photographs since 1976 (Fig. 5), and the length change record at Franz Josef Glacier8 (Figs 2 and 5).

Overall, 1,035 direct glacier mass balance measurements made at Franz Josef and Tasman glaciers between 1972 and 2011 are used to tune and evaluate the glacier mass balance model. In all, 455 of these mass balance measurements on Franz Josef Glacier between 1 April 2000 and 31 March 2002 are used in model tuning, while 580 measurements from before, and after, this tuning period provides an independent assessment of model performance. Overall, agreement between modelled and measured mass balance is excellent (Supplementary Fig. 2 and Supplementary Table 2). The comparison of modelled and measured accumulation also provides a secondary test of the precipitation input data.

Additional model runs were conducted to ensure the results from our diagnostic experiments were not unduly influenced by model parameter choices; the range of output values that resulted from the range of parameter choices given in Supplementary Table 1 is the blue and red uncertainty bands for precipitation and temperature, respectively, in Fig. 6d. Repeat model simulations using statistically downscaled re-analysis temperature (rather than interpolated station data) showed similar patterns (Supplementary Fig. 4) and fall within the uncertainty bounds generated using different model parameter sets (Fig. 6d).

To evaluate the influence of our simple debris-cover ablation-reduction scheme and our assumption of glacier surfaces always being at the melting point, we tested this simple scheme against our full subsurface and debris conduction calculations using the same model experiments (Supplementary Fig. 4). We found that the more complete thermal calculations gave a very similar result to our simplified analysis. Including conductive melt under debris, and subsurface heat fluxes in snow and ice, the contribution of temperature variations to volume changes is generally decreased, while the contribution of precipitation variations to volume changes is close to zero over the full period. The contribution of other variables also changes, but remains minor.

Climate analysis

We analysed interannual ice volume changes in the ‘standard run’ in order to link this signal to changes in the atmosphere and oceans. First, the control run, which contains the overall trend of volume loss over the study period, based on the mean climate, is used to de-trend the mass balance series to remove the effect of the large imbalance between present-day climate and glacier extent (shown by the generally negative mass balance in Fig. 4a). The year-to-year glacier mass balance volume anomaly variations (Fig. 5a) are ranked from largest to smallest values. The upper and lower quintiles from this ranking are then used to identify the strongest 8 years for ice volume gain and loss (Supplementary Table 4).

We calculated means and variances for several climate indices that correspond to these individual glacier years. These glacier years are further subdivided into four seasons; April, May, June (AMJ), July, August, September (JAS), October, November, December (OND) and January, February, March (JFM). The first half of the glacier year (i.e., the cool season; April to September), and the second half of the glacier year (i.e., the warm season; October to March) are also considered. Two-sample t-tests are applied to determine whether the means for the two groups are significantly different at the 95% confidence level (Supplementary Figs 6 and 7). The climate indices considered in our analysis include the Southern Oscillation Index, the Niño 3.4 index, New Zealand Trenberth Indices, regional blocking indices for the southern mid-latitudes, the Southern Annular Mode, the Pacific South American mode and Zonal Wave 3 (ref. 35).

We used the Australian Bureau of Meteorology data set ( http://www.bom.gov.au/climate/current/soihtm1.shtml) for the Southern Oscillation Index73, while the Niño 3.4 (N3.4) index was sourced from a National Oceanic and Atmospheric Administration data set ( http://www.cpc.ncep.noaa.gov/data/indices/ersst3b.nino.mth.ascii). The Trenberth M1 index is the normalized surface pressure difference between Hobart, Australia and Christchurch, New Zealand, and is a broad measure of meridional flow over the South Island and Tasman Sea51. The Trenberth MZ3 index is the normalized surface pressure difference between New Plymouth (North Island, New Zealand) and Chatham Island to the east of New Zealand and is a broad measure of South West/North East flow over the South Island of New Zealand51. Southern Annular Mode, Zonal Wave 3 and Pacific South American pattern index anomaly values were calculated relative to the 1981–2010 climatological normal interval using climate re-analysis data30. Both Zonal Wave 3 and the Pacific South American pattern are based on monthly mean 500-hPa geopotential height; the Pacific South American index is the amplitude time series of the second Empirical Orthogonal Function calculated over the southern extratropics south of 20°S, and Zonal Wave 3 is calculated according to ref. 35. The Southern Annular Mode index is calculated using normalized monthly zonal mean sea-level pressure difference between 40°S and 65°S according to ref. 74. The Southwest and Southeast Pacific Blocking Indices were calculated following ref. 75. We also analysed the monthly frequencies of occurrence of 12 New Zealand region ‘Kidson’ synoptic types76 in relation to glacier variability.

Seasonal composite spatial patterns for sea surface temperature, temperature at the 850 hPa level, geopotential height anomaly at 500 and 1,000 hPa, zonal wind anomalies and precipitable water content anomalies are constructed for April to June, July to September, October to December, January to March, cool (April to September) and warm (October to March) seasons, and glacier years (April to following March) using online re-analysis from the National Oceanic and Atmospheric Administration ( http://www.esrl.noaa.gov/psd/cgi-bin/data/composites/printpage.pl). All of the re-analysis data are available for the period of our mass balance simulation (1972–2011). Anomalies relative to 1981–2010 are calculated for all climate variables (Figs 7 and 8 and Supplementary Fig. 5). Sea ice concentration data extend back only to late 1981, and we utilized the National Oceanic and Atmospheric Administration monthly mean ice concentration data set available at ( http://www.esrl.noaa.gov/psd/cgibin/db_search/DBSearch.pl?) to generate spatial maps of sea ice concentration anomalies. Sea ice anomalies in the individual seasons in the composite plots are relative to the 1982–2010 average, and two seasons in the gain quintile (1975–1976 and 1979–1980) are therefore omitted from the compositing analysis because of temporal limitations of that data set. In addition, a sea ice concentration index for the Ross Sea sector of the Southern Ocean are utilized for cluster analysis, available from the National Snow and Ice Data Center ( http://nsidc.org/data/smmr_ssmi_ancillary/area_extent.html).

Evaluating natural and anthropogenic climate influences

A recent global-scale attribution study of glacier mass balance changes3 concluded that New Zealand glacier changes were consistent with mass balance model simulations forced by global climate models that included all forcings (natural and anthropogenic), rather than mass balance model simulations forced by natural forcings alone. However, a limitation of the application of this study to New Zealand glaciers3 is the use of a global data set of glacier mass balance observations to assess their models53. These ‘observations’53 are 5-year averages, interpolated from global measurements during periods when no local mass balance observations are available in New Zealand (for example, the 1980s and 1990s). In this paper we provide a New Zealand-specific indication of mass balance variations for these unmeasured periods.

We compare the full ensemble and ensemble means of simulations that included all forcings, and just natural forcings, in this global attribution study3 to those of our ‘standard run’ mass balance simulation to assess whether our simulations also reveal this anthropogenic influence on New Zealand glacier mass balance (Fig. 9 and Supplementary Table 6). We used two-sample t-tests to assess whether the means of these mass balance simulations were significantly different. We also used consecutive t-tests on 11-year arrays, starting from 1973 to 1983 and ending at 2001–2011 to assess whether this relationship changed through time.

Data availability

Additional data sets not included in the Supplementary Information are available from the corresponding author on reasonable request.