A geomorphology based reconstruction of ice volume distribution at the Last Glacial Maximum across the Southern Alps of New Zealand

We present a 3D reconstruction of ice thickness distribution across the New Zealand Southern Alps at the Last Glacial Maximum (LGM, c. 30 e 18 ka). To achieve this, we used a perfect plasticity model which could easily be applied to other regions, hereafter termed REVOLTA (Reconstruction of Volume and Topography Automation). REVOLTA is driven by a Digital Elevation Model (DEM), which was modi ﬁ ed to best represent LGM bed topography. Speci ﬁ cally, we removed contemporary ice, integrated offshore bathymetry and removed contemporary lakes. A review of valley in- ﬁ ll sediments, uplift and denudation was also undertaken. Down-valley ice extents were constrained to an updated geo-database of LGM ice limits, whilst the model was tuned to best- ﬁ t known vertical limits from geomorphological and geochronological dating studies. We estimate a total LGM ice volume of 6,800km 3 , characterised pre- dominantly by valley style glaciation but with an ice cap across Fiordland. With a contemporary ice volume of approximately 50 km 3 , this represents a loss of 99.25% since the LGM. Using the newly created ice surface, equilibrium line altitudes (ELAs) for each glacier were reconstructed, revealing an average ELA depression of approximately 950m from present. Analysis of the spatial variation of glacier-speci ﬁ c ELAs and their depression relative to today shows that whilst an east-west ELA gradient existed during the LGM it was less pronounced than at present. The reduced ELA gradient is attributed to an overall weakening of westerlies, a conclusion consistent with those derived from the latest independent climate models. © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
The New Zealand Southern Alps is a key site for glacial reconstruction studies and is one of only three mid-latitude southern hemisphere locations at which extensive Pleistocene glaciation occurred. Furthermore, the geomorphological record of the Last Glacial Maximum (LGM) in New Zealand has excellent spatial and temporal coverage  and has long been recognized for its potential for improving understanding of former climatic conditions via glacier geometry reconstruction (e.g. Porter, 1975b).
LGM ice thickness has previously been reconstructed for individual catchments in New Zealand (e.g. Putnam et al., 2013) and at the regional scale (Golledge et al., 2012), although uncertainty still exists regarding the extent, timing and associated climatic conditions.
Glacial reconstruction in New Zealand presents significant challenges as the extreme relief and high precipitation of parts of the Southern Alps gives rise to a complex yet large-scale system (Chinn et al., 2014). To date, the work by Golledge et al. (2012) is the only attempt at a regional reconstruction of LGM ice thickness, and that was achieved using the Parallel Ice Sheet Model (PISM). They used a Digital Elevation Model (DEM) and climate parameter inputs (spatially uniform scaling of contemporary conditions) to best fit the areal ice limits delimited by Barrell (2011) (Fig. 1). Some relatively large (tens of kilometers) discrepancies must be noted however, where the PISM was unable to recreate known downvalley limits derived from field-based observations and independently dated landforms. Specifically, in comparison to the mapping of Barrell (2011), Golledge et al. (2012) underestimated ice extent in the eastern outlets including the Rakaia and Waimakariri (Fig. 1).
As an alternative to top-down climate models, a number of bottom-up, 'topography based' techniques have been developed for estimating contemporary ice thickness (Linsbauer et al., 2012;James and Carrivick, 2016) and palaeo-ice thickness (Pellitero et al., 2016). With the Southern Alps presenting excellent input datasets for LGM ice thickness reconstruction in this manner (high resolution DEMs and well constrained down-valley ice limits), and with climate driven models currently unable to resolve known ice limits in some catchments (i.e. Golledge et al., 2012), we developed and applied a topography driven model to produce a new estimate of distributed LGM ice thickness. Our approach is similar to that described by Pellitero et al. (2016) and in effect a reverse of the VOLTA contemporary ice thickness model developed by James and Carrivick (2016). As such the model is hereafter referred to as REVOLTA (Reconstruction of Volume and Topography Automation).
Using this approach, the aim of this study is to reconstruct the Fig. 1. Location map of regions referred to in this study and previous delineations of LGM ice extent. Note the mismatch between extents estimated by Barrell (2011) and Golledge et al. (2012). Inset: Modelled ice thickness and extent by Golledge et al. (2012) and extent mapped by Barrell (2011) for the Rakaia Valley.
LGM distributed ice thickness across the Southern Alps of New Zealand to help improve our knowledge of glaciological and associated climatic conditions during the LGM in the Southern Hemisphere. Being far from the centres of Northern Hemisphere glaciation, research in New Zealand provides the potential to examine the climatic effects of changes in Northern Hemisphere ice sheets on Southern Hemisphere climate, as well as events of origin in the Southern Hemisphere .
2. Methods: a framework for reconstructing distributed ice thickness REVOLTA initially calculates the surface profiles of former glaciers along a network of centrelines using an exact solution of a 'perfect plasticity' approach. Ice thickness values are iteratively calculated along each centreline, starting at the former terminus position. A 'nearest neighbour' routine is used in conjunction with the DEM to convert the centreline thickness points into a 3D distributed thickness. REVOLTA expands on the 2D only model presented by Benn and Hulton (2010) and is similar in nature to Pellitero et al. (2016). A schematic diagram of the REVOLTA framework is shown in Fig. 2 and major processing steps are described in the following sections.

Constructing a DEM to represent the LGM
A digital elevation model (DEM) forms the primary input of REVOLTA. Whilst contemporary DEMs are often used to represent former glacier beds (Trommelen and Ross, 2010), the most accurate results will be achieved if the DEM is representative of the period of reconstruction (i.e. the LGM). To achieve this, we used a contemporary 8m resolution DEM (Geographix, 2012), to which we applied a number of corrections: (i) removing contemporary ice thickness, (ii) integrating offshore bathymetry, (iii) removing modern-day lakes that fill the remnants of troughs evacuated by LGM ice. The impact of uplift, denudation and postglacial sediment valley infill was also considered.

Removal of contemporary ice from the DEM
The Southern Alps currently supports around 3,000 glaciers covering an area of approximately 1,200 km 2 (Chinn, 2001). Contemporary ice was removed from the DEM to estimate a surface which better represents LGM bed conditions. This removal was achieved using the distributed ice thickness results of the VOLTA model  to compute modern day ice thickness (e.g. as for Southern South America by Carrivick et al. (2016) and as for the entire Antarctic Peninsula by ). VOLTA was run on all contemporary glaciers >1 km 2 and resulted in a total (contemporary) ice volume of 40.68 km 3 with a maximum thickness of over 550 m. Full details of the VOLTA model can be found in James and Carrivick (2016), whilst details of the contemporary ice thickness distribution of the Southern Alps can be found in James (2016).

Appending an offshore DEM
Sea level during the LGM was approximately 125 m below current levels (Milne et al., 2005), with palaeoglaciers on the West Coast of New Zealand inferred to extend into these currently submarine environments. Specifically, analysis of areal LGM limits delineated by Barrell et al. (2011) revealed that 2.5% of a total LGM glaciated area occurred beyond the present coastline and consequently outside the limits of the standard DEM. Therefore, an alternative topographical source was required for these presentlyoffshore areas. For this purpose, we used the National Institute of Water and Atmospheric Research (NIWA) 250 m resolution bathymetric dataset (Mitchell et al., 2012) and these data were appended to the existing DEM. To create a seamless dataset, the bathymetry was resampled to the resolution of the land surface (8 m) and merged.

Contemporary lakes
Whilst palaeoglaciers extended past the contemporary coastline, they also occupied areas where contemporary lakes now exist, for which standard DEMs represent the water surface level and not the underlying bathymetry, thus obscuring the former glacier bed topography. Substantial lakes inboard of LGM moraine sequences of the Southern Alps are over 450 m deep (e.g. Irwin, 1980), highlighting the importance of considering their bathymetry. Many of these lakes developed at the end of the LGM and were in contact with a calving ice margin (Sutherland et al., 2019).
To produce a DEM more representative of LGM conditions, the water depth of major lakes within the model domain was subtracted from the initial DEM. Each map was georeferenced, with water depth contours ranging in interval between 5 m and 100 m digitized. Bathymetry was adjusted to the lake level of the DEM before a gridded surface was interpolated using the hydrologically correct ANUDEM algorithm (Hutchinson, 1989). These data were subsequently merged with the DEM, helping to produce a surface more representative of LGM bed conditions.
Whilst sediment deposition has occurred in these lakes since the LGM, there is insufficient data to remove this from the DEM. Of the limited data which does exist, cores from Lake Ohau found sediment thicknesses (deposited since lake formation at the end of the LGM) of 82 m and 43 m beneath 101 m and 68 m of water respectively (Levy et al., 2018). This shows that whilst the presence of such sediment needs to be acknowledged, water depth alone accounts for the majority of the difference between the LGM bed and the contemporary DEM surface. By removing the water depth, our study accounts for over 60% of the difference between the LGM bed and the contemporary DEM, something which has not been considered by previous studies (e.g. Golledge et al., 2012).

Uplift and denudation
Uplift and denudation are processes that modify the landscape, meaning their influence on the DEM needs to be considered. Whilst these changes are commonly too small to introduce significant error (Finlayson, 2013) and are thus commonly not included in ice thickness modelling studies (Golledge et al., 2012;McKinnon et al., 2012), here we assess how robust their omission is in terms of impacting ice volume of the Southern Alps.
Post-LGM isostatic rebound in New Zealand is considered to have been minimal (Pickrill et al., 1992), with Mathews (1967) suggesting a maximum of 30 m rebound. However, tectonic uplift rates in the Southern Alps are high, with the Pacific plate being compressed into the Australian plate (Coates, 2002). Whilst long term Cenozoic uplift rates of 0.1e0.3 mm yr À1 were estimated by fission track thermochronology (Tippett and Kamp, 1995), rates during the Pleistocene and Holocene are thought to have been greater (Adams, 1980). The Paringa River site (location shown in Fig. 1) provides one of the few robust constraints, with Pleistocene to Holocene rates calculated at 8 ± 1 mm yr À1 (Norris and Cooper, 2001). However, this rate is not thought to be representative of the Southern Alps as a whole as the uplift pattern was highly spatially variable, as evidenced by tilted lake shorelines (Adams, 1980).
Since AD 2000, GPS measurements have been collected along a transect running broadly perpendicular to the Alpine Fault, providing further information on the distribution of uplift rates across the Southern Alps. Results show a maximum average vertical rate of 5.4 mm yr À1 for the transect (Beavan et al., 2007). However, and critically, rates reduce rapidly with distance from the Alpine Fault, with measurements 29 km from the fault recording an uplift rate of just 1.8 ± 1 mm yr À1 , whilst no significant vertical rise was measured at 68 km (Beavan et al., 2007). As these GPS measurements record current rates, it is acknowledged that they are not necessarily applicable to longer (glacial) timescales. Nonetheless, with contemporary uplift rates reducing exponentially with distance from the Alpine Fault and no significant rise recorded at 68 km from the fault, and in lieu of any other geological data pertaining to long-term and spatially-distributed uplift rates, it is likely that much of the model domain for an LGM ice reconstruction would fall within relatively low uplift areas. This interpretation is consistent with the inferred uplift pattern of Adams (1980), with the majority of the LGM ice extent occupying low uplift regions.
Importantly, denudation effectively counteracts uplift, with long term rates in the Southern Alps approximately equal (Adams, 1980). For example, whilst the highest peak in the Southern Alps (Aoraki, Mount Cook) is 3,724 m high, the total amount of uplift during orogeny is estimated to be approximately 15e20 km (Kamp et al., 1989). Tippett and Kamp (1995) provide one of the few quantifications of denudation, using fission track thermochronology to estimate long term Cenozoic denudation rates ranging from 2.5e0.5 mm yr À1 , decreasing with distance from the Alpine Fault. This is supplemented by contemporary river load measurements by Adams (1980) who estimated the total river load of the Southern Alps to be 700 ± 200 Â 10 9 kg yr À1 whilst basin averaged rates for the Nelson/Tasman region are estimated to be between 112 and 298 t km À2 yr À1 (Burdis, 2014). Whilst providing a useful insight into denudation rates, these estimates are not sufficient to construct a regional distribution required for DEM modification.
In summary, whilst it is acknowledged that uplift and denudation have modified the geomorphological expression of the landscape, the approximate counterbalance between the two processes limits their impact between the land surface at the LGM and the contemporary DEM. Furthermore, with rates of both uplift and denudation decreasing with distance from the Alpine Fault, much of our model domain falls within relatively low uplift/denudation areas. With the New Zealand land mass largely in its present configuration by the mid-Pleistocene (Barrell, 2011), we therefore made no further modifications to the DEM. These issues of tectonic deformation and denudation will become increasingly prevalent for older glaciations, so care must be taken if seeking to reconstruct earlier glacial episodes.

Postglacial sediment valley infill
Large volumes of sediment accumulated in valley floors is a commonly noted feature of the Southern Alps (McKinnon et al., 2012). This process of landscape modification has previously been recognized as a potential source of uncertainty when reconstructing palaeoglaciers in New Zealand (Golledge et al., 2012) and needs to be considered when using a DEM to represent a former glacier beds.
It was initially postulated that the majority, or entirety of valley fill in the New Zealand Southern Alps was deposited post-LGM (Suggate, 1965;Adams, 1980) simply due to the 'absence of any conflicting evidence ' Adams (1980, p. 77), with some contemporary research supporting the hypothesis that glaciers did indeed erode to bedrock at the LGM (Thomas, 2018). However, this notion has been challenged for the down valley reaches of LGM glaciers in the Southern Alps, with Rother et al. (2010) suggesting that recent glaciations in the Hope Valley (see Fig. 1 for location) were not powerful enough to cause extensive erosion and overrode pre-LGM sediments rather than removing them. At this location, infrared stimulated luminescence (IRSL) dating and sediment analysis found the majority of material was deposited during the pre-LGM period of 95.7 to 32.1 ka with only approximately 30 m of sediment accumulation occurring post LGM, overriding 200 m of previous deposits . This pattern is consistent with seismic investigation in the Franz Josef Valley by Alexander et al. (2014), suggesting that the LGM bed was well above the bedrock. A modelling study by McKinnon et al. (2012) also supports the concept of LGM glaciers overriding pre-existing sediment, using a variety of models to simulate the LGM bed of the Pukaki glacier. The notion of pre-LGM sediment surviving in glaciated valleys is also consistent with studies outside of New Zealand, with pre LGM deposits found beneath more recent sediments in the European Alps (e.g. Hinderer, 2001).
In summary, recent field based observations and modelling studies provide evidence that LGM glaciers often overrode preexisting sediments to some extent and did not always reach bedrock. Whilst it may seem surprising that gravels were able to survive the entirety of the LGM, Shulmeister et al. (2019) suggests that rapid water drainage through these stratified gravels resulted in motion being initially limited to the slow and inefficient process of internal ice deformation, allowing ice advance without substantial erosion of the underlying substrate. As such, whilst techniques exist to estimate and remove the entirety of sediment infill from a DEM (e.g. Harbor and Wheeler, 1992;Jaboyedoff and Derron, 2005), this approach is not deemed appropriate for the regional scale modelling of LGM glaciers in the New Zealand Southern Alps.
With some sediment surviving, LGM bed elevations lie somewhere between the contemporary DEM and the level that would exist if all postglacial valley-fill were removed. However, since inferred LGM sediment thickness is only available for a few select locations (e.g. McKinnon et al., 2012), there are not enough data to generate a regional distribution of the LGM bed surface beneath Holocene valley fill sediments. Therefore, with the consensus in the literature suggesting that at least some sediment was present pre-LGM, for the regional scale modelling, REVOLTA was applied without removing post LGM sediment.

Distributed post-LGM sediment thickness of the Pukaki Valley.
As discussed in section 2.1.5, regional modelling was performed without removal of any post-LGM sediment as there is insufficient data to generate a regional distribution. However, modelling by McKinnon et al. (2012) provides sufficient information to estimate post-LGM sediment thickness for the Pukaki Valley and consequently allowed us to test the sensitivity of REVOLTA at this location. Specifically, McKinnon et al. (2012) found that a mass-flux balance model (as described by Anderson et al. (2004)) produced the most reliable LGM bed profile and associated post-LGM sediment thicknesses. The reader is directed to McKinnon et al. (2012) and Anderson et al. (2004) for a full discussion of the model choice and its mechanics.
Spatially distributed post-LGM sediment thickness was estimated by firstly identifying the areal extent of sediment in-fill and subsequently interpolating the sediment thickness in this region using the transect data of McKinnon et al. (2012). The areal extent was identified using a modified version of the algorithm presented by Straumann and Korup (2009), using a GIS derived stream network as 'seed cells' and iteratively assessing the relative elevation of neighboring cells. Specifically, we set a gradient threshold of 1.5 between neighboring cells as this was found to best match the visual geomorphological evidence. The mathematics of the algorithm can be found in Straumann and Korup (2009) and for our purposes we coded the process in the Python programming language. Once the areal extent of sediment infill was delineated, we interpolated the post LGM thickness within this region using the mass flux balance bed profile generated by McKinnon et al. (2012). Interpolation was achieved using the ANUDEM algorithm as this is designed specifically for topographical applications (Hutchinson, 1989).

Areal extent of LGM glaciers
Alongside the DEM, areal ice extents are required as a model input for REVOLTA, with the downstream point of each centreline defining the starting point for each glacier reconstruction. The latest revision of LGM ice extent mapping was produced by Barrell (2011), including data based on the QMAP (Quarter Million Scale Mapping) project, a nationwide geological mapping project where the ages of quaternary deposits were mapped in terms of marine isotope stages (Heron, 2014). As such, we use the mapping of Barrell (2011) as the input for REVOLTA, with some minor adjustments made to account for new and unexploited evidence, as summarised in Table 1. The final extent outline we used can be seen in Fig. 1 and a full discussion for each individual adjustment can be found in James (2016).

Centreline generation and point ice thickness estimation
REVOLTA initially estimates ice thickness values at regular intervals along a network of centrelines in a similar manner to other 2D (Benn and Hulton, 2010) and 3D (Pellitero et al., 2016) approaches. Centrelines for all major catchments and tributaries were constructed using a GIS based hydrological routing approach, as used in other glaciological studies (e.g. Schiefer et al., 2008). Centrelines were manually checked and adjusted to best define the centre of each glacier branch.
Ice thickness values were subsequently calculated at 50 m intervals along each centreline, using a perfect plasticity approach similar to that developed by Benn and Hulton (2010) and used by Pellitero et al. (2016). The original basis for this method was introduced by Schilling and Hollin (1981), stating that: Where h i is the ice surface elevation, H is the ice thickness, T b is the basal shear stress (equal to the yield stress), f is a shape factor representing the proportion of the driving stress supported by the bed, p is ice density (917 kg m À3 ), g is gravitational acceleration (9.81 m s À2 ) and i is a specified interval of distance along the centreline (50 m in this study).

Shear stress
Basal shear stress (T b ) is a critical parameter of Equation (1) as it determines at which point deformation will occur. T b is generally between 50 and 150 kPa for mountain valley glaciers (Paterson, 1994), although it may vary between individual glaciers due to various factors (e.g. basal sliding, ice viscosity, subglacial deformation). A constant shear stress value was used along the length of each centreline as this has been shown to adequately reproduce ice thickness estimates along the length of a centreline (Li et al., 2012) and has been successfully used in other palaeoglacier reconstruction studies (Rea and Evans, 2007). We calibrated T b using independent estimates of ice thickness from a variety of published records. This method is more reliable than using globally averaged values (e.g. Locke, 1995;Rea and Evans, 2007) and is the preferred choice for calibrating perfect plasticity based models (Benn and Hulton, 2010;Pellitero et al., 2016).
Although the standard perfect plasticity model assumes motion by internal deformation only, basal sliding may also contribute to motion if glaciers are wet based. It is unclear to what extent basal sliding occurred during the LGM in New Zealand, with Golledge et al. (2012) suggesting the geomorphological evidence noted by Mager and Fitzsimons (2007) (striated and abraded bedrock, deformed subglacial till and the transport of erratics through subglacial pathways) as evidence of its occurrence. Conversely, Shulmeister et al. (2019) propose that water rapidly drained through the stratified gravel beds, minimizing basal sliding. In either case, basal sliding is accounted for in this study as T b is calibrated from local landform evidence, effectively prescribing a softer ice rheology, as described by Veen (2013).
Whilst prescribing a softer ice-rheology in this manner is suitable for reconstructing former profiles, it should be noted that T b effectively becomes a 'lumped parameter' and further calculations, such as estimates of ice velocity may be spurious. Basal sliding also needs to be carefully considered if T b is calculated differently (e.g. from globally averaged values).

Shape factor
It is unrealistic to assume that the bed shear stress equals the driving stress in mountain glaciers as friction from the valley walls provides resistance to flow (Paterson, 1994). This effect is commonly incorporated into perfect plasticity models using a 'shape factor' (f ) as shown in Equation (1). REVOLTA dynamically adjusts f depending on the local valley geometry using a modified version of the method described by Benn and Hulton (2010): Where A is half the glacierised area (of the cross section), p is half the glacierised perimeter and H is the ice thickness at the centreline. Initially, an approximate ice surface elevation and cross section is generated by using a standard f value of 0.8 (Nye, 1965). From this initial cross section, a 'valley geometry based' f value can be calculated using Equation (2) which is used to update the ice surface elevation. In a similar approach to that of Benn and Hulton (2010), an iterative process is applied to determine if the valley geometry based f value adequately describes the updated cross section (defined as within an error of ±0.1), recalculating the ice surface elevation and valley geometry based f value until an appropriate solution is found.

Creation of an ice surface
For each point at which the ice surface elevation is estimated along the centreline, a 'nearest neighbour' approach is used to find the closest point on each of the valley sides which is of equal or greater elevation. The procedure is executed twice; once for each valley side, and the resultant nearest points effectively mark the lateral extent of the glacier being reconstructed. These points on the valley sides are assigned a zero ice thickness and are used with the centreline point ice thickness values to generate a fully 3D distributed ice surface.

Automated palaeo-ELA estimation
The LGM ice surface generated by REVOLTA was used to estimate Equilibrium Line Altitudes (ELAs) during the LGM. Firstly, individual glacier catchments were approximated by generating watersheds from the newly created LGM ice surface. ELAs for each individual glacier were then calculated using the Area-Altitude Balance Ratio (AABR) method (Osmaston, 2005). This was achieved using Python code modified from Pellitero et al. (2015) and as used by  for Little Ice Age glacier ablation areas, but in this case using the newly generated LGM ice surface and outlines as inputs. With insufficient data (i.e. mass balance) available to parameterize the balance ratio directly, we used representative balance ratios suggested by Rea (2009), as recommended by Pellitero et al. (2015). To assess the sensitivity of the specific balance ratios used, we calculated the ELAs under two 10 km upstream (to Mt. Iron moraine) IRSL dating (Wyshnytzky, 2013) and radiocarbon dating (McKellar, 1960). Milford Sound 3 km upstream (to correspond with submarine terminal moraine) Mapping and SEDs (Dykstra, 2012).
Waitutu region Adapted to extent of Ward (1988) Mapping (Ward, 1988) scenarios; firstly using a 'global' value of 1.75 for all glaciers and secondly using a 'mid-latitude maritime' value (1.9) for glaciers to the west of the main divide and a continental mid-latitude value (1.59, derived from the European Alps) for those to the east. The values in the latter scenario were chosen to best accommodate the east-west climate gradient of the Southern Alps (Henderson and Thompson, 1999), which is also believed to have been a prominent during the LGM (Drost et al., 2007;Shulmeister et al., 2019).

Contemporary lake bathymetry
Table 2 details the lakes removed from the DEM and summary statistics. The location of each lake can be seen on Fig. 1. The lakes have a combined volume of 240.6 km 3 , highlighting the importance of their consideration for palaeo ice thickness estimation. Furthermore, this is the first known attempt to quantify the volume of these lakes, revealing that Lake Wakatipu is volumetrically the largest (64.2 km 3 ). This is interesting as the title of the 'largest lake of the South Island' is sometimes noted as Lake Te Anau as it is the largest by surface area (e.g. McLintock, 1966).

Shear stress
The vertical extent of LGM glaciation in New Zealand is undefined for the majority of individual catchments. The available evidence, where geomorphological mapping (e.g. Barrell et al., 2011) has been combined with geochronological dating and some numerical modelling, is summarised in Table 3. This dataset permitted multiple model runs of varying shear stress to be made of each LGM glacier to manually best-fit modelled palaeo ice surfaces to that field-measured and dated. Any surface exposure dates cited hereafter have been recalibrated to the updated production rate from Putnam et al. (2010) in a similar manner to that of Darvill et al. (2016).
We aimed to generate a profile slightly above the dated landform evidence due to the general convex profile of glaciers in the ablation zone (Nesje, 1992) where centreline ice elevations are above those on the valley sides, with a differential of 50 m approximated in previous research deemed appropriate (Glasser et al., 2011).
The results of the shear stress parameterisation are shown for each location in Fig. 3. A value of 30 kPa best-fitted the available evidence, providing an acceptable fit to all locations. Since these catchments include examples from both the east and west of the Main Divide, it appears that for the purposes of our study, a single regionally-averaged value of 30 kPa is appropriate for the entire Southern Alps. Whilst it is acknowledged this approach does not accommodate some inter-catchment variability, we consider it appropriate for our modelling because our aim is to provide geometric and volumetric details, rather than information on local ice dynamics or thermal regime, for example. Fig. 4a shows the results of the reconstructed LGM ice thickness across the Southern Alps, with an estimated total volume of 6,800 km 3 . The REVOLTA model produces large, predominantly valley-constrained, glaciers flowing to the east of the Main Divide (e.g. Rakaia and Pukaki Glaciers, Fig. 4b and c) and some piedmont style glaciers coalescing on the West Coast (e.g. Haast region, Fig. 4d). With the LGM shoreline estimated 125m below present (Milne et al., 2005), only 1 glacier extends marginally past the LGM coastline (Fig. 4d), accounting for just 0.5 km 3 of ice (<0.007% of the total modelled LGM ice volume).

Distributed ice thickness and total ice volume
A large number of nunataks are present, especially around the high peaks of the central region. Far fewer nunataks are present south of the Lake Wakatipu region and indeed REVOLTA produces a localised ice cap in the Fiordland region (Fig. 4e). A maximum modelled ice thickness of 1,108 m is located at the Wakatipu Glacier (see Fig. 4 for location).

Sensitivity analysis e shear stress
We conducted a sensitivity analysis of T b using the LGM Pukaki Glacier as a test case. Analysis was performed for the Pukaki Glacier due to the abundance of dated landform evidence (Table 3). As previously noted, the 'best fit' T b value of 30 kPa was purposely chosen to generate a profile slightly above the dated landform evidence. This is due to the general convex profile of glaciers in the ablation zone (Nesje, 1992) where centreline ice elevations are above those on the valley sides, with a differential of 50 m approximated in previous research deemed appropriate (Glasser et al., 2011). Considering such a convex profile, a 'low' T b estimate (15 kPa) was used to approximate the lowest elevation landform evidence whilst a 'high' T b estimate (40 kPa) was used to exceed the elevation of the geomorphological evidence by over 50 m. Fig. 5 shows the impact of altering T b within the bounds described: a T b of 15 kPa (Fig. 5a) produces an area 17.5% less than the 'best fit' 30 kPa scenario whilst the 40 kPa scenario (Fig. 5c) produces an area 9.5% greater than the baseline. For ice volume there was an underestimate of 47.4% (15 kPa model) and an overestimate of 30.5% (40 kPa model). Much of these differences are driven by the specific topography of the Pukaki Valley, although they provide a first order estimate of the model sensitivity.

Sensitivity analysis e post LGM sediment infill
Sensitivity analysis was also conducted with regards to post LGM sediment infill. Analysis was performed for the Pukaki Valley as this is the only location of the Southern Alps where sufficient data on sediment thickness exists (McKinnon et al., 2012). Using the methods described in Section 2.1.5.1, Fig. 6a shows the estimated post-LGM sediment thickness distribution (and areal extent) within the Pukaki Valley, revealing up to 348 m of post-LGM infill. Modelled LGM ice thickness is shown using the original (no sediment removed) DEM in Fig. 6b and for the newly created 'sediment removed' version in Fig. 6b. Corresponding cross sections of the LGM bed profile, original DEM and ice surfaces are shown in Fig. 6d. For this modelling, all other parameters were set as per the standard model (T b of 30 kPa). When considering overall ice volume in the study region, modelling with the standard DEM predicts a volume 4.0 km 3 (7.3%) lower than that with the sediment removed.

ELA reconstruction
To compare LGM ELAs with contemporary ELAs, we use a digitized version of the contemporary ELA trend surface produced by Chinn and Whitehouse (1980). We use this surface as our baseline as it is effectively an ELA 0 surface as most glaciers were close to zero balance during this period (Chinn et al., 2005). Table 4 summarises the results of the reconstructed LGM ELAs under both balance ratio scenarios (as described in section 2.7) alongside the contemporary values for corresponding regions, whilst Fig. 7 graphically displays the results for the balance ratio 1.9 west/1.59 east scenario. LGM ELAs were substantially lower than present, with regionally averaged estimates of approximately 900 m. There was large spatial variability of ELAs during the LGM, as there is at present (Carrivick and Chase, 2011), with LGM glaciers flowing to the west coast exhibiting ELAs approximately 320 m lower than those to the east. A gradual south west e north east trend can also be observed (Fig. 7a), with the lowest values found in the south west. Superimposed upon these trends is further local variability, with localised ELA low points found where additional mountain masses are parallel and offset to the west of the main divide (e.g. Solution Range, Fig. 7a) whilst higher ELAs are found where low passes penetrate the main divide (e.g. Hollyford Pass, Fig. 7a). Both balance ratio scenarios produce similar results, although a marginally stronger east-west differential and associated ELA gradient are evident under the 1.9 west/1.59 east scenario.
LGM ELAs were depressed by approximately 950 m on average, although depressions were greater for glaciers to the west of the main divide (under both balance ratio scenarios). Calculation of average ELA gradients across three transects (as originally used by Chinn and Whitehouse (1980)) reveals that east-west ELA gradients during the LGM were markedly reduced at locations where they are presently very high (Fig. 7, T2: Mt. Cook and T3: Olivine Range transects) whilst increasing at the northernmost transect (T1: Whitcombe Pass) where a lower ELA gradient is currently present.

Style and nature of glaciation during the LGM
Our estimate of 6,800 km 3 ice volume for the Southern Alps during the LGM is similar to that by Golledge et al. (2012) of 6,400 km 3 (adjusted herein for the same spatial extent). The comparison is especially interesting because of the markedly different approaches: Golledge et al. (2012) took climate proxy data to drive an ice flow dynamics model, whereas we utilised geomorphological evidence of glacier extent to derive distributed glacier ice thickness independent of climate estimates. With current ice volume of approximately 50 km 3 (Chinn, 2001), the Southern Alps has lost approximately 99.25% of its glacial ice since the LGM, highlighting the climatic sensitivity of the glaciers. For comparison, Greenland has lost approximately 36% of its ice since the LGM and Antarctica just 7% (Bintanja et al., 2002).
There is some debate and uncertainty surrounding the vertical ice extent in New Zealand during the LGM although there is a growing consensus that glaciers were somewhat constrained by topography and a full ice cap was not present. Anderson and Mackintosh (2006) suggest the abundance of extremely sharp ridgelines as evidence they were not overridden by ice whilst modelling of individual catchments clearly shows exposed ridgelines around the Mt. Cook region (Putnam et al., 2013). Although the previous regional modelling by Golledge et al. (2012) presents a continuous 'icefield' style of glaciation, ice thickness over many of the peaks and ridgelines is actually modelled to be extremely low.
The results of REVOLTA in this study are consistent with a valley constrained style of glaciation, with exposed ridges and nunataks, especially in central and northern regions (Fig. 4). The main Table 3 Locations used for manually parameterising shear stress value for use in REVOLTA. See Fig. 1  exception to this is in the Fiordland region, where REVOLTA suggests a localised ice cap style of glaciation (Fig. 4e). Interestingly, this pattern is almost identical to that proposed by Chinn et al. (2014) who postulated that the Darran Range (adjacent to Milford Sound) marked an approximate boundary between a southern ice cap and a northern region where peaks and ridges penetrated the LGM ice.
The REVOLTA modelling approach also provides an indirect insight into potential basal conditions as the parameterisation scheme for T b effectively accounts for basal sliding (Veen, 2013). Specifically, we found a T b value of 30 kPa best represented landform evidence, with a range between 15 kPa and 40 kPa deemed plausible. If motion was solely by internal deformation we would expect T b to be around 50e150 kPa (Nye, 1952;Paterson, 1994). As such, we interpret that basal sliding was prominent for NZ glaciers during the LGM. Whilst this is in agreement with other studies (e.g. Golledge et al., 2012) due to sedimentary and geomorphological evidence (Mager and Fitzsimons, 2007;Evans et al., 2013), Shulmeister et al. (2019) suggests this is not always the case. In their study, it is suggested that basal substrate plays a crucial role, with rapid water drainage through stratified gravels limiting motion to internal deformation. Conversely, when the bed consists of a high mud content (e.g. lake basins), rapid basal sliding may occur. Whilst our study supports the notion of generally wet and sliding conditions at the regional scale, further research is needed to quantify variability over space and time.
Whilst the total volume and most of the longitudinal profiles presented in Fig. 3 are similar to those of Golledge et al. (2012), the Milford profile (Fig. 3e) is markedly different, with differentials of ice surface elevation in excess of 300 m in places. The good agreement between REVOLTA and the (albeit limited) surface exposure dates in this catchment suggests that the results of REVOLTA are the most plausible in this location. One potential reason for Golledge et al. (2012) overestimating glacier ice surface elevations and volume in this specific catchment (compared the age controlled geomorphological evidence and REVOLTA) is their assumption of spatially-uniform differentials of temperature and precipitation between the LGM and present. This is noted in their analysis, where scenarios which better reproduced known glacial limits to the east of the main divide (e.g. the Rakaia) produced overextended glaciers in the central west coast areas. In essence, their 'optimal reconstruction' can be seen as a compromise between under-estimating certain glaciers to the east of the main divide (as discussed in their analysis and shown in Fig. 1) and over-estimating those of the central west coast (as shown in Fig. 3e).

ELA reconstruction and implications for climatic conditions
Our LGM ELA estimates for individual glaciers of between 375 m and 1,502 m are of a similar range to those reported by Golledge et al. (2012) (approximately 350 me1,450 m) and show an overall similar pattern. The two balance ratio scenarios we used produced an almost identical regionally averaged ELA (900 m compared to 898 m) whilst the 1.9 west/1.59 east scenario produced a slightly increased east/west differential, as would be expected (Table 4). As the differences are minimal, we suggest that either scenario is appropriate for the purposes of our modelling, although the 1.9 west/1.59 east scenario may be preferred as this better represents the known conditions of the glaciers currently in the region (Chinn et al., 2005). Importantly, the global scenario (balance ratio for all glaciers ¼ 1.75) still reproduces a strong west-east differential, so we can be confident that this differential is not due to the choice of balance ratios on each side of the main divide.
A high east-west ELA gradient has long been recognized between contemporary Southern Alps glaciers and is often interpreted as a response to the current precipitation pattern (Porter, 1975a). As such we infer the presence of an east-west ELA gradient during the LGM as evidence that westerly circulation was in operation to some extent during the LGM, a view shared by other authors (e.g. Rother and Shulmeister, 2006). Crucially however, our   . 7) shows that LGM ELA gradients were markedly reduced compared to present (reducing from 31 to 14 m/km À1 and 20 and 12.5 m/km À1 respectively) whilst the currently low ELA gradient at the Whitcombe Pass (3.7 m/km À1 ) was enhanced during the LGM to 12.6 m/km À1 . Our assertion of spatial variability in climatic change is in line with a growing body of evidence from other studies using a variety of techniques, although to our knowledge this study is the first to use glacial reconstruction for such a purpose in New Zealand. McKinnon et al. (2012) collated estimates of temperature deviation (for the New Zealand LGM) from different studies, finding a large range in reported values. For example, Sikes et al. (2002) found that sea surface temperature south of Chatham rise were 8 C cooler than present whilst Marra et al. (2006) inferred very moderate cooling in the Canterbury region, as little as 2 C on average (albeit probably during a short inter-stadial within the extended LGM). The notion of spatial variability is supported by climate modelling by Drost et al. (2007) who estimated that temperature reduction varied spatially between 2.5 C and 4 C from pre-industrial levels.
Although LGM precipitation estimates are less well constrained than temperature, there is evidence to suggest that both the amount and distribution were significantly different from present. Palaeo-vegetation studies suggest a major reorganization of vegetation patterns during the New Zealand LGM, with a general consensus of drier conditions overall (McGlone et al., 1993). Precipitation patterns were also likely to be different during the LGM, with recent research by Shulmeister et al. (2019) postulating changes driven by a shift in the track of westerlies. Precipitation modelling by Drost et al. (2007) also supports the concept of a generally drier climate, albeit with substantial local variation.
Most recently, downscaled climate modelling by Karger et al. (2017) and Brown et al. (2018) has produced high resolution (~1 km) gridded estimates of LGM mean annual temperature and precipitation, which were downloaded from www.paleoclim.org. Fig. 8a displays the mean annual temperature distribution at the LGM using the CCSM4 model, whilst Fig. 8b shows the corresponding values at present (Landcare Research Ltd, 2003b). Similar maps are shown for precipitation distribution, with Fig. 8c showing the modelled LGM distribution and Fig. 8d showing the contemporary distribution (Landcare Research Ltd, 2003a). By differencing the LGM estimates and corresponding contemporary datasets, the differentials can be calculated for both temperature (Fig. 8e) and precipitation (Fig. 8f). Note that LGM glaciers are delineated in Fig. 8e as no correction has been made for changes in ice surface elevation and associated effect of lapse rate. Importantly, the downscaling precipitation algorithm used by Karger et al. (2017) and Brown et al. (2018) incorporates orographic predictors including wind fields, valley exposition, and boundary layer height, with a subsequent bias correction. We are therefore confident that the modelled precipitation data (Fig. 8c) adequately captures the magnitude of orographic precipitation enhancement in the highrelief setting of the Southern Alps.
These new, high resolution datasets form an independent source for assessing the spatial variability of climatic parameters. Fig. 8e and f demonstrates that changes in temperature and precipitation distribution varied spatially, as predicted by the ELA reconstructions of this study (Fig. 7). Fig. 8e shows the majority of the Southern Alps experienced cooling in the range of 2 e 5 C, broadly similar to the New Zealand average of 4.6 C reported by Drost et al. (2007) although perhaps slightly lower than the 6.5 C suggested by Golledge et al. (2012). Whilst the negligible cooling estimated in north-west regions is surprising, this is within the ranges reported in the literature, with Marra et al. (2006) suggesting a temperature increase of up to 0.5 C during summer LGM interstadials at Lyndon Stream (see Fig. 8e for location) whilst Samson et al. (2005) reported winter sea surface temperature decreases of just 0.9 C (offshore of New Zealand North Island).
For precipitation, Fig. 8f shows major changes in the amount and distribution during the LGM. Whilst an east-west precipitation gradient is clearly still modelled during the LGM (Fig. 8c), this is much reduced, with less precipitation on the west coast, which currently has some of the highest rainfall on Earth (Henderson and Thompson, 1999). The pattern is consistent with the assertion of an overall drier conditions during the LGM (McGlone et al., 1993;Drost et al., 2007). Drost et al. (2007) also found the greatest precipitation decreases in Western regions, with average reductions of 30% (from pre-industrial levels) in Southern Otago-Southland regions.
In essence, the climate model estimates based on the work of Karger et al. (2017) and Brown et al. (2018) are in correspondence with the results from the ELA reconstructions in this study (Fig. 7a): namely a weakened (albeit still important) east-west climatic Fig. 7. a) LGM ELAs (AABR, balance ratio 1.9 west/1.59 east). b) LGM ELA depression from present. gradient with substantial local level variability. The agreement is especially interesting because of the markedly different approaches: Karger et al. (2017) and Brown et al. (2018) used a global climate model to estimate temperature and precipitation wheras we utilised geomorphological evidence to derive ice thickness and associated ELAs.

Conclusions
In this study we used geomorphological and geochronological records from the Southern Alps of New Zealand to drive a glacial reconstruction model (REVOLTA) for the LGM. We modelled a total volume of 6,800 km 3 of glacier ice across the Southern Alps during the LGM, representing a 99.25% volume loss to present. Analysis of the distributed ice thickness reveals a 'valley-constrained' style of glaciation with many exposed ridges and nunataks, even when parameterised to fit the maximum vertical extent of geomorphological evidence. An ice-cap style of glaciation was predicted for the Fiordland region, where maximum surface elevations are lower. We found a low shear stress value of approximately 30 kPa was required to best fit the geomorphological evidence, from which we infer that at least some glaciers were warm-based during the LGM.
Reconstructed ELAs reveal a regionally averaged depression of approximately 950 m during the LGM. This ELA depression is consistent with findings of other studies. However, and crucially, we found substantial spatial variability in ELA depression, with east-west ELA gradients during the LGM generally decreased. In conjunction with evidence from the latest independent climate models (Karger et al., 2017;Brown et al., 2018), we propose that westerly circulation was decreased during the LGM compared to at present, resulting in decreased precipitation to the west of the main divide.
The recognition of spatially variable climatic differentials is important because they have not been previously considered when reconstructing LGM ice thickness in New Zealand. Namely, the only previous such model by Golledge et al. (2012) assumed spatiallyuniform differentials of temperature and precipitation between the LGM and present. Whilst producing an acceptable fit at the regional scale, they were unable to resolve known ice limits in some catchments (by up to 40 km). We suggest that some of these anomalies may be due to the use of unrealistic uniform climate differentials, especially precipitation distribution to which glaciers of the Southern Alps are known to be highly sensitive (Rowan et al., 2014). The advent of independent high resolution estimates of palaeo climatic conditions (Karger et al., 2017;Brown et al., 2018) may help to address some of these issues, allowing further differences between 'top down' climate based reconstructions and 'bottom up' landform based reconstructions (e.g. REVOLTA) to be assessed in more detail. The framework presented here could easily be applied to other formerly glaciated regions, helping to improve our knowledge of the global climate system. https://doi.org/10.1016/j.quascirev.2019.06.035.