Early glacier advance in New Zealand during the Antarctic Cold Reversal

Glacial landscapes preserve records of past climate change. Investigating the glacier–climate system over the Late Quaternary provides information about past climate change and context for present‐day glacier response to climate warming. Using 28 beryllium‐10 (10Be) surface exposure dates and snowline reconstructions, we present glacier fluctuations and climate changes for the Antarctic Cold Reversal in the Ahuriri River catchment, Southern Alps of New Zealand (44°7′50″S, 169°38′29″E). Prominent terminal and lateral moraine features from the upper right tributary of the Ahuriri River valley have exposure ages of 14.5 ± 0.3, 13.6 ± 0.3 and 12.6 ± 0.2 ka, suggesting retreat of the glacier during the Antarctic Cold Reversal. Maximum elevation of lateral moraines (MELM) and accumulation area ratio (AAR) suggest snowline elevations at these ages were ≤700, ≤630 and ~360 m lower than today, respectively. This equates to air temperatures ≤3.9, ≤3.5 and 2.3 ± 0.7 °C lower than today (1981–2010), assuming no changes in past precipitation. Ice‐sculpted bedrock surfaces bound by a lateral moraine at nearby Canyon Creek have an age of 13.1 ± 0.3 ka, indicating the moraine correlates with those in the Ahuriri upper right tributary. MELM and AAR reconstructions from the Canyon Creek suggest that snowline elevations at 14.5–13.6 ka were ≤500 or ~380 m lower than today, corresponding to air temperatures ≤2.8 or 2.4 ± 0.7 °C lower than the present‐day (1981–2010). Our results provide insight into the structure of the Antarctic Cold Reversal in the Southern Alps, showing that the largest glacier advance occurred at the start of this interval at c. 14.5 ± 0.3 ka and was followed by gradual retreat. We hypothesize that the early cooling and glacier readvance in New Zealand at the onset of the Antarctic Cold Reversal were triggered by a latitudinal shift of the Southern Hemisphere westerly wind belt.


Introduction
Worldwide glacial retreat during the last major deglaciation coincided with a large-scale reorganization of oceanic and atmospheric circulation patterns, global sea-level rise and changes in greenhouse gas concentrations (Mix et al., 2001). This period of time, known as the Last Glacial Termination (or last deglaciation) extends from the end of the Last Glacial Maximum (LGM, 19-20 ka;Clark et al., 2009) to the onset of the Holocene (11.7 ka; Denton et al., 2010). The later stages of the last deglaciation between 15.0 and 11.5 ka are known as the Lateglacial interval. During the Lateglacial, the warming trend recorded in Antarctic ice cores was interrupted by the Antarctic Cold Reversal (ACR, 14.5-12.5 ka) (Jouzel et al., 1995), which was subsequently also recognized in the southern mid-latitudes, including New Zealand (e.g. Putnam et al., 2010a;Kaplan et al., 2010;Kaplan et al., 2013;Pedro et al., 2016). Climate fluctuations during the Lateglacial are important because they provide an example of abrupt, nonlinear climate change during a past period of global warming, which may have relevance for present-day climate change .
Mountain glaciers are highly responsive to climate fluctuations and provide crucial information on regional and global climate (Oerlemans, 2005). Cosmogenic exposure dating techniques using nuclides such as 10 Be are a widely used tool to understand past glacial fluctuations (Balco, 2011(Balco, , 2020. Glacier snowlines are also frequently used to extract quantitative information about paleo-temperature and precipitation for former glaciers . Glacier records from New Zealand situated in the southwestern Pacific provide a rare record of Southern Hemisphere climate changes, and are ideal tools for testing hypotheses about the competing drivers of Lateglacial climate change (Kaplan et al., 2010).
In the Southern Alps, where the present work focuses, only a small number of Lateglacial cosmogenic moraine chronologies have been documented. These records provide evidence that a glacier-advance event during the Lateglacial culminated at 13 ka (e.g. Kaplan et al., 2010;Putnam et al., 2010a) (see Fig. 1 for site locations), suggesting close agreement with Antarctic ice core temperature records that indicate warming abruptly resumed at this time (Jouzel et al., 2001;Pedro et al., 2011). Some of these moraine records also include older ages that suggest glaciers may have been slightly larger at around 14.5 ka  and 14.1 ka (Putnam et al., 2010a). However, due to the limited number of ages available, the significance of these data for understanding millennial-scale climate events of the Lateglacial remains uncertain. Further records of glacier change from wellpreserved moraine sequences are required to provide more complete information about the structure of mid-latitude climate variability during this millennial-scale event.
Here, we target glacial landforms of Lateglacial age at two different sites in the Ahuriri River valley in the Southern Alps of New Zealand. These two sites were chosen because moraines and glaciated bedrock sequences are well preserved. Furthermore, the former glaciers occupied single valleys of simple geometry, avoiding potential complications associated with tributary glaciers, making these basins ideal for reconstructing past snowlines (e.g. Pellitero et al., 2015).

Southern Alps
The Southern Alps of New Zealand extends~700 km from northeast to the southwest, forming an orographic barrier to the prevailing westerly winds in the southwestern Pacific region. Mean annual precipitation ranges from 3000 to 10 000 mm on the western slopes of the Alps, declining to about 1000 mm on the eastern slopes of the mountain range (Henderson and Thompson, 1999). Mean annual air temperatures range from 9 to 13°C at sea level across the South Island. Temperatures decrease at a rate of about~1°C for every 200 m increase in altitude (Norton, 1985). There are 17 peaks that exceed 3000 m in height, and the range culminates in its central section at Aoraki/Mount Cook (3724 m a.s.l.). Over 3000 glaciers cover~800 km 2 of the Southern Alps (Baumann et al., 2020). Snowline altitude ranges from 1500 m a.s.l. just west of the Main Divide at Franz Josef Glacier to >2200 m a.s.l. east of the Main Divide at the Ben Ohau Range (Lamont et al., 1999).
The Ahuriri River originates on the eastern slopes of the central Southern Alps, and forms the border between the Canterbury and Otago regions (44°16′56″S, 169°36′32″E). While few precipitation measurements exist for the Ahuriri River valley, it probably receives high precipitation in the upper parts of its catchment near the Main Divide, and lower precipitation in its lower parts, which extend eastwards.
The general geographical setting of the Ahuriri River valley was described by Tielidze et al. (2021Tielidze et al. ( , 2022. In the current study we focus on two different sites in the valley, the upper right tributary of the Ahuriri, and Canyon Creek (Fig. 2).

Upper right tributary of the Ahuriri River valley
A small tributary valley with a length of~2.7 km connects to the upper western side of the Ahuriri River valley (44°7′41″S, 169°37′52″E). The main part of this valley is surrounded by steep, eroded slopes (especially the northern bench) resulting in extensive debris deposits (colluvium) on the lower valley walls and floor. The upper section of the valley reaches its maximum altitude at~2360 m a.s.l., while the lowest part descends to an elevation of~1150 m a.s.l. The lowest portion of the valley is represented by a glacial riegel (rock bar) that has been exposed by ice erosion. The surface of this riegel is fractured and partly vegetated.
Several lateral and terminal moraines exist at different altitudes in the lower part of the valley. Two clearly visible lateral moraine ridges hereafter called outer-1 (44°7′57″S, 169°38′34″E) and outer-2 (44°7′55″S, 169°38′30″E) run parallel to each other,~20 m apart, on the southern flank of the valley (Fig. 2). These moraines project into the main valley, indicating that the former glacier exited the  (Sime et al., 2013). White arrows indicate direction of the SWW belt. The yellow, blue and white dotted lines represent positions of the Polar Front, Sub-Antarctic Front and Sub-Tropical Front, respectively (Darvill et al., 2016). The red dots are locations of key ice core palaeoclimate records illustrated in Figs. 11 and 13 and mentioned in the text [Dronning Maud Land (DML), Byrd, Siple Dome, Talos Dome, Dome C, Law Dome] (Jouzel et al., 2001;Bereiter et al., 2015;Pedro et al., 2011;Stenni et al., 2011). (b) Location map of South Island of New Zealand showing study sites mentioned in the text and the Ahuriri River valley, the focus of this study. See Fig. 2 for a detailed view of the selected field sites. [Color figure can be viewed at wileyonlinelibrary.com] tributary and terminated in the Ahuriri River valley at the time of moraine formation. An additional well-preserved terminal moraine ridge (hereafter called inner) is situated 140 m lower in elevation relative to the outer ridges 1 and 2 representing a thinner and shorter glacier that terminated just inboard of the riegel (44°7′54″S, 169°38′42″E). There are further ridges of uncertain origin running parallel to each other on the northern bench of the lower valley slope (44°7′44″S, 169°38′52″E) that were not investigated during the field campaign.

Canyon Creek
Canyon Creek is one of the largest tributaries of the Ahuriri River. The upper western section of this valley is surrounded by Huxley Range with mountain peaks reaching~2200 m a.s.l.
The Huxley Range also separates Canyon Creek from Hunter River (Lake Hāwea basin) to the west. Other mountain peaks at the crest of the western surrounding ridge reach~1800 m a.s.l. The largest extant glacier of the Ahuriri River catchment (Thurneysen Glacier -0.9 km 2 ) is located in the headwater of this valley (44°9′56″S, 169°35′55″E). Separately, a large cirque exists on the western headwater slope of the Canyon Creek valley. The floor of this cirque contains a small lake (tarn) bound by fresh moraines that probably indicate late Holocene deglaciation. Several small alluvial fans and debris cones are observed mainly on the right side (west relative to Canyon Creek) slopes across the river.
A large ice-sculpted bedrock outcrop is situated on the floor of the middle section of Canyon Creek and cut through by the river. The largest exposure is represented on the true left side (44°11′22″S, 169°35′28″E) (Fig. 2). Several small ponds occur on the surface of this bedrock while the rest of the area is partly vegetated. The elevation of this bedrock ranges from 1250 to 1500 m a.s.l.
A prominent linear left-lateral moraine ridge is situated on the eastern slope of the valley (44°11′28″S, 169°35′34″E) just above the mentioned bedrock outcrop (Fig. 2). This ridge is almost 800 m long and ranges in elevation from 1305 to 1500 m a.s.l. The orientation of the moraine ridge is not parallel to the valley (or to the river bed), but the crest of the moraine projects towards the valley centre, suggesting that the palaeo glacier terminated not far down-valley from this moraine. A smaller (~180 m length) and less visible fragment of a right-lateral moraine ridge is preserved on the opposite slope, ranging from 1230 to 1275 m a.s.l. (44°11′30″S, 169°34′47″E). These paired lateral moraines are separated by glacial riegel of~200 m relief.

Field work
We targeted ice-sculpted bedrock surfaces (in both the upper right tributary, and Canyon Creek) and greywacke boulders on the well-preserved moraine ridges (only from the upper right tributary) for cosmogenic surface exposure dating. Many of the boulders contained quartz veins of various thicknesses (5-30 mm). Overall, 28 samples were collected using a hammer and a chisel during the austral summer of 2020 (Fig. 3). We preferentially sampled horizontal rock surfaces to ensure optimal exposure of the sampling sites to the incoming cosmic-ray flux. The thickness of the extracted samples ranged from 1.2 to 7 cm (Table 1). We used a geological compass and clinometer to measure the azimuthal elevations of the surrounding landscape. These were combined with strike and dip observations of each boulder surface to calculate topographic shielding of the cosmic ray flux using the online exposure (formerly known as CRONUS-Earth) calculator (Balco et al., 2008). A Trimble GeoXH differential GPS system was used to measure the position of each boulder including altitudes. GPS data were differentially corrected by post-processing using the Land Information New Zealand PositioNZ network to achieve sub-metre accuracy in each dimension. Sample elevations are reported relative to the 2016 New Zealand Geoid. Boulders were measured, described and photographed. All field observations are given in Table 1.

Exposure age calculations
Age calculations were carried out using the online exposure (formerly known CRONUS-Earth) calculator, version 3 (Balco, 2011; available: https://hess.ess.washington.edu/). We assumed a rock density of 2.7 g/cm 3 and applied sample thickness corrections based on calliper measurements made prior to crushing (Table 1). We use the 'Macaulay' calibration data set of Putnam et al. (2010b) due to its close geographical proximity to our samples. For comparison, we also present ages using the primary global production rate calibration dataset of Borchers et al. (2016). We present ages using the St (Stone, 2000), Lm (Balco et al., 2008) and LSD (Lifton et al., 2014) scaling schemes. However, we only discuss ages based on the 'Lm' scaling scheme. This decision was made to provide consistency between this and previous studies from the Southern Alps which used the same 'Lm' scaling scheme (e.g. Putnam et al., 2010a;Kaplan et al., 2013;Tielidze et al., 2022), as well as the similarity of results across these scaling schemes for our sample sites. In addition, the chi-squared outlier detection routine was used from the exposure age calculator version 3 (Balco, 2017a(Balco, , 2017b to assess if spread in the 10 Be exposure ages for a single moraine belt is consistent with them having been deposited synchronously with each other. Where sample sets for a single landform fail the chi-squared test for a single population, we follow the iterative procedure described by Balco (2017a) to identify and remove anomalous exposure ages. Both study sites lie above the elevation (~1200 m) at which snow lies persistently during the winter months in the Southern Alps, raising the possibility that exposure to cosmic radiation at our samples may have been partially shielded by seasonal snow cover. Kaplan et al. (2010) observed that snow cover at~1500 m a.s.l. in Irishman Stream basin (~35 km from our study sites) has rarely been thicker than 2 m or so over the last few decades allowing salient boulders to stand out prominently from the snow blanket. Most of the boulder heights from our study site range from~1 to 3 m and they are located on the crest of the moraine ridges, and thus we are confident that snow shielding at these sites is negligible relative to other uncertainties (e.g. measurement, production rate scaling). For bedrock samples, snow cover may be more significant, yet the magnitude and duration of seasonal cover remains unconstrained. Theory indicates that a 1 m of average density snow (0.3 g/cm 3 ) covering a rock surface for 4 months of the year reduces 10 Be production via spallation by~5% (Gosse and Phillips, 2001), which equates to 500-600 years for our bedrock sample sites. We present ages without this speculative correction, but acknowledge that the apparent ages may underestimate the true age by a few centuries. Overall, bedrock samples comprise only a small portion of our data and our assumptions outlined above do not alter the main conclusions of our study.

Reconstruction of glacier geometry, snowline and temperature
Accurate reconstruction of past equilibrium-line altitudes (ELAs) (commonly called 'snowlines') require accurately determining the extent and morphology of former glaciers. ELA reconstructions, along with the determined age of the reconstructed glaciers, provide a quantitative proxy for palaeoclimatic conditions (e.g. Porter, 1975;Benn et al., 2005). Detailed geomorphological mapping (e.g. Tielidze et al., 2021) is required for reconstruction of past glacier extent, including terminal moraines and trimlines, which allow the shape of former glaciers to be most accurately reconstructed (Benn et al., 2005). Based on the geomorphological map by Tielidze et al. (2021) and high-resolution Google Earth imagery, we reconstructed the surface of the former glaciers in Canon Creek and the upper right tributary of the Ahuriri. Maximum upstream limits of former glaciers are defined by the catchment head, while the downstream and marginal ice limits are defined precisely by terminal moraines where they are preserved. For the palaeo glacier area uncertainty we used a buffer method which is broadly adopted for modern glacier mapping (e.g. Tielidze and Wheate, 2018). A buffer width of 50 m was created along the glacier outline, and the uncertainty term was calculated as an average ratio between the original glacier area and the area with a buffer increment. This generated an average uncertainty of the mapped glaciers area of ±9%. Contour lines derived every 50 m enabled us to create 30-m-resolution digital elevation models (DEMs). The area between each pair of successive contours was measured automatically in ArcGIS by the ELA calculation toolbox (Pellitero et al., 2015). An accumulation area ratio (AAR) assumes that the accumulation area of a glacier occupies a fixed proportion of the total glacier area (Benn and Gemmell, 1997). Knowledge of mass-balance gradients is not required for this method, but a reconstructed 2D glacier surface is needed to calculate the AAR. An AAR value of 2:1 or 0.67/0.33 is frequently used for snowline estimation in New Zealand (e.g. Porter, 1975;Kaplan et al., 2010;Eaves et al., 2016Eaves et al., , 2017, and was also recently used for LGM snowline reconstruction in the Ahuriri River valley (Tielidze et al., 2022). Here we implement this technique using the toolbox developed for ArcGIS by Pellitero et al. (2015) with an AAR ratio of 0.67 ± 0.05 (1 SD).
We also used the Maximum Elevation of Lateral Moraines (MELM) method as a supplementary tool for snowline estimation. This was due to the absence of a terminal moraine that could potentially accompany the prominent lateral moraine ridges in Canyon Creek (44°11′28″S, 169°35′ 34″E) and in the upper right tributary of the Ahuriri River valley (outer-1, 44°7′57″S, 169°38′34″E and outer-2, 44°7′ 55″S, 169°38′30″E). The MELM method is based on the assumption that the highest lateral moraine ridges constrain the minimum altitude of a former snowline, because net annual glacier ablation, and hence ice-marginal deposition, only occurs below the snowline. Thus, the former snowline must lie above the highest elevation of the lateral moraines (Benn and Lehmkuhl, 2000). However, sometimes it is difficult to assess whether moraine deposition started immediately down-glacier of the snowline or whether a lateral moraine is preserved entirely in the upper part. Thus, again, only a lower limit for the snowline can be estimated (Atle, 1992). The MELM method does not require accurate topographic information for the whole glacier but only accurate determination of the highest elevation of the lateral moraines. In this study, the highest point of a lateral moraine feature was taken as the lower limit of palaeo snowline in the Canyon Creek and in the upper right tributary of the Ahuriri River valley (outer-2).
For consistency with our previous study from the Ahuriri River valley (Tielidze et al., 2022), we use a present-day snowline elevation of~2000 m a.s.l. for Thurneyson Glacier (44°9′55″S 169°35′50″E). This value closely matches with a recently reported snowline elevation of 1957 m a.s.l. for Thurneyson Glacier based on four-decades (1981-2010) of snowline observations from the National Institute of Water and Atmospheric Research (Lorrey et al., 2022). Using the present-day temperature lapse rate this allowed us to reconstruct the corresponding temperatures from Lateglacial time. The present-day temperature lapse rate in New Zealand centres on~5°C km −1 (Norton, 1985;Tait and Macara, 2014), although this is poorly constrained by observations. Thus, we follow Eaves et al. (2016) in selecting temperature lapse rates from an evenly distributed range of −4 to −7°C km −1 . The algorithm repeatedly resamples these input parameters providing a population of pELA and dT estimates that represent the full window of parameter uncertainty. This approach was also successfully used by Tielidze et al. (2022) for constraining temperature anomalies in the Ahuriri River valley during the LGM and termination. The approach that we apply to reconstruct past temperatures assumes that temperature changes are the main driver of glacier changes in the Southern Alps, and specifically that precipitation changes and other meteorological variables (e.g. cloud cover) play a minor role. While this assumption is a simplification, it is justifiable given the evidence from mass balance sensitivity studies on present-day glaciers in the Southern Alps (e.g. Anderson et al., 2010;Anderson and Mackintosh, 2012) as well as sensitivity analyses conducted in palaeoglaciological modelling studies that target a similar time period Eaves et al., 2017). We revisit this assumption and its influence on our findings in the Discussion.

Results
The spatial distribution of moraine and bedrock features across the upper right tributary of the Ahuriri River valley and Canyon Creek suggests the occurrence of several glacial re-advance or stillstand periods during the Lateglacial. The timing of these glacial phases is constrained by 28 10 Be surface exposure ages spanning from 14.5 ± 0.3 to 12.6 ± 0.2 ka (here and henceforth all quoted uncertainties are 1-sigma). Sample details for surface exposure ages from all moraine systems are given in Table 1, while cosmogenic 10 Be exposure ages are listed in Table 2 and Figs. 4-6. All age calculations are referenced to in calendar years before sample collection (2020). The exact data blocks that were entered into the online calculator for age calculation are provided in Supporting Information Tables S1-S5.

The upper right tributary of the Ahuriri River valley
Outer moraine ridge-1 Outer moraine-1 is a single ridge,~170 m long at an elevation ranging from 1230 to 1300 m a.s.l. Four exposure ages from this moraine ridge are tightly clustered and range from 14.2 ± 0.4 to 14.8 ± 0.3 ka (Figs. 4 and 5). These four ages yield an internal error-weighted mean age of 14.5 ± 0.3 ka (n = 4; no outliers).

Outer moraine ridge-2
Outer moraine-2 is inboard of outermost moraine-1. It is much longer than outer-1 stretching to about 560 m. The crest of this moraine feature ranges in altitude from~1160 to~1370 m a.s.l. and is scattered by large 1-to 3-m (diameter) sized greywacke boulders. Six samples from this lateral moraine feature returned 10 Be ages from 13.3 ± 0.3 to 15.2 ± 0.5 ka. The ages from outer moraine ridge-2 were grouped around an exposure age of 13.6 ± 0.3 ka (n = 5; 1 outlier: AU-26-76), based on the internal error-weighted mean of the six exposure ages (Table 2 and Figs. 4 and 5).

Inner moraine ridge
The crest of the inner terminal moraine feature ranges from 1160 to 1230 m a.s.l. Eleven samples from this moraine ridge show individual apparent exposure ages ranging from 2.0 ± 0.2 to 13.0 ± 0.3 ka (Table 2 and Fig. 4). Two boulders (AU-27-91 and AU-27-95) with the youngest ages of 2.0 ± 0.2 and 11.0 ± 0.3 ka are not situated on the main moraine ridge. We consider that these boulders fell from a nearby slope and came to rest on an older deposit. Thus, we do not include these ages in our calculation of the mean landform age. The remainder of the nine samples from this terminal moraine ridge show individual apparent exposure ages ranging from 12.3 ± 0.3 to 13.0 ± 0.3 ka. The ages from this moraine landform were grouped around an exposure age of 12.6 ± 0.2 ka (n = 9; two outliers), based on the internal error-weighted mean of the 11 exposure ages.

Riegel-bedrock near the Inner moraine ridge
A glacial riegel (rock bar) is situated at the lowest portion of the valley at an elevation of 1150 m a.s.l. One bedrock sample (AU-26-90b) from this riegel was collected about 80 m downvalley from the crest of the inner moraine ridge. This sample with an age of 13.8 ± 0.5 ka is stratigraphically consistent with the age of the inner moraine ridge and thus provides some constraint on the timing of glacier retreat and exposure of the riegel. The age is consistent with the abandonment of outer moraine ridge-2 (13.6 ± 0.3 ka), suggesting that the glacier retreated from outer moraine ridge-2 to a position inside the tributary valley at this time.

Canyon Creek
Six exposure ages from the bedrock surface situated on the left side (east relative to Canyon Creek) of the river (44°11′22″S, 169°35′28″E) range from 10.5 ± 0.2 to 13.6 ± 0.3 ka. The ages from this bedrock surface have a mean exposure age of 13.1 ± 0.3 ka (n = 6; one outlier -CC-27-98b), based on the internal error-weighted mean of the six exposure ages (Table 2 and Fig. 5). The Canyon Creek bedrock ages suggest that this site was covered by ice during lateral moraine formation (outer-1, 14.5 ± 0.3 ka; outer-2, 13.6 ± 0.3 ka) at the upper right tributary of the Ahuriri River valley.
Palaeo glacier geometry, snowline and temperature estimation

Upper right tributary of the Ahuriri River valley
We reconstructed the geometry of the former glacier when the terminus was at the inner latero-terminal moraine position at 12.6 ± 0.2 ka in the upper right tributary of the Ahuriri River valley. Based on this reconstruction, the glacier had a total surface area of 1.8 ± 0.2 km 2 while its length was~2.5 km. The  (1 and 2). These observations indicate that the snowline was ≤700 m lower than present-day snowline at 14.5 ± 0.3 ka and ≤630 m lower than present-day snowline during 13.6 ± 0.3 ka, suggesting that local air temperature was ≤3.9 and ≤3.5°C lower than present respectively.

Canyon Creek
Sampled bedrock surfaces from Canyon Creek are situated immediately inboard of a prominent left lateral moraine ridge (44°11′28″S, 169°35′34″E) which must have been deposited before the bedrock became exposed at 13.1 ± 0.3 ka (Fig. 6). Thus, we tentatively suggest that this moraine probably correlates with the outer-1 and outer-2 moraine ridges from the upper right tributary of the Ahuriri River valley, which were deposited at 14.5-13.6 ka. Consequently, these moraine features from both Canyon Creek and the upper right tributary are likely to indicate the timing of formation (re-advance or stillstand) of both palaeo glaciers (Canyon Creek and upper right tributary) at 14.5-13.6 ka.
Based on this geomorphic and age relationships between those two sampled sites we reconstructed the geometry of the palaeo glacier from Canyon Creek. During the interval 14.5-13.6 ka the glacier had a total surface area of 6.0 ± 0.8 km 2 while its length and thickness was~5.5 km and~300 m respectively. The estimated snowline based on the AAR method was~1620 ± 60 m a.s.l., while the lowest limit of palaeo snowline based on the MELM method was at about ≥1500 m a.s.l. These correspond with a snowline depression of~380 m (AAR) or ≤500 m (MELM) relative to the presentday, which represents a temperature of −2.4 ± 0.7 or ≤2.8°C lower than present during the moraine formation in the Canyon Creek at 14.5-13.6 ka (Figs. 9 and 10), assuming no change in precipitation.

Moraine chronology
Our surface-exposure chronology from the upper right tributary moraines of the Ahuriri River basin indicates that an advance of the palaeo glacier culminated at 14.5 ± 0.3 ka (outer moraine ridge-1), while the next re-advance or stillstand phase occurred at 13.6 ± 0.3 ka (outer moraine ridge-2). About 1 kyr later (12.6 ± 0.2 ka), the former glacier had lowered ~140 m and built another prominent terminal moraine ridge in the lower section of the valley. We did not estimate geometry of the palaeo glacier at 14.5 ± 0.3 and 13.6 ± 0.3 ka from this valley since both lateral moraine features (1 and 2) do not indicate the terminal position of the former glacier, but only the lateral-marginal locations. The 10 Be constraints from the Canyon Creek valley indicate the prominent, valley-floor bedrock outcrop became free from ice cover at 13.1 ± 0.3 ka. This result is consistent with the upper right tributary where our ages also recorded ice thinning and glacier retreat at this time (i.e. between outer-2 and inner moraines).

Early glacier advance in the Lateglacial
Previous studies in the Southern Alps, for example at Irishman Stream (Kaplan et al., 2010) and Birch Hill (Putnam et al., 2010a), argue that the coldest part of the Lateglacial was at 13.0 ka, at the end of the ACR (Fig. 11d,e). However, a few studies indicate that the maximum Lateglacial advance or stillstand was earlier. For example, Kaplan et al. (2013) provide ages for a glacier advance at 14.5 ± 0.4 ka (n = 4) from two adjacent moraine ridges (two ages per moraine ridge) (Fig. 11e). Two boulders of similar age have also been constrained at the Birch Hill outboard moraine ridge (14.1 ± 0.4 ka, n = 2) ( Fig. 11d)  .
Our data clearly demonstrate the structure of this early Lateglacial event in the Southern Alps, with a large glacier advance evident at the start of this interval at c. 14.5 ± 0.3 ka, a small decrease in ice extent between 14.5 ± 0.3 and 13.6 ± 0.3 ka, a larger recession between 13.6 and 12.6 ka, and one more re-advance or stillstand at 12.6 ka (Fig. 11g). Consequently, our study provides confirmation that a glacier re-advance or stillstand occurred at 14.5 ± 0.3 ka, suggesting that the coldest part of the Lateglacial reversal was earlier than previously thought in New Zealand (e.g. Kaplan et al., 2010;, assuming no changes in precipitation. Our finding is comparable with that from outlet glaciers in Torres del Paine (51°S, south Patagonia, Chile) which reached their maximum extent at 14.2 ± 0.6 ka (García et al., 2012).
Our results from outer moraine ridge-2 (13.6 ± 0.3 ka) and inner terminal moraine ridge (12.6 ± 0.2 ka) show consistent long-term trends with existing dated glacier landforms and snowline reconstructions from the Southern Hemisphere. For instance, our 10 Be constraints from the outer moraine ridge-2 (13.6 ± 0.3 ka) are consistent with 10 Be ages of 13.3 ± 0.3 ka at Mid-Macaulay moraines (middle ridge, Tekapo basin) in the Southern Alps (Putnam et al., 2010a). Our ages are also in good agreement with a glacier record from Monte San Lorenzo (47°S), central Patagonia, which show that the Río Tranquilo Glacier re-advanced at 13.9 ± 0.7 and 13.7 ± 0.5 ka (Sagredo et al., 2018). The youngest Lateglacial advances in the upper right tributary of the Ahuriri River valley (inner terminal moraine -12.6 ± 0.2 ka) overlap in age with the 12.7 ± 0.5 ka at the intermediate moraine in the outer part of Irishman Stream (Ōhau basin) in the Southern Alps (Kaplan et al., 2010). This is also consistent with a moraine deposit dated at 12.4 ± 0.3 ka from Lacteo valley in Monte San Lorenzo, central Patagonia (Mendelová et al., 2020).

Early cooling in the Lateglacial
Our geologically constrained estimates of ELA changes, relative to the present, afford quantitative insights of climate changes associated with past glacier geometries. Glacier mass balance is primarily determined by temperature and precipitation, while other variables such as solar radiation, wind speed and relative humidity also contribute to varying degrees depending on the climatic setting (Oerlemans, 2001;Mackintosh et al., 2017). In New Zealand, previous studies have shown that air temperature changes dominate over other Figure 11. (a) δD (deuterium) and temperature anomalies (Jouzel et al., 2001;Röthlisberger et al., 2002) from Dome C, Antarctica. (b) Composite CO 2 record from Antarctica (Bereiter et al., 2015). (c) Biogenic opal flux (wind-driven upwelling) (Anderson et al., 2009). (d,e) Normal kernel density plots from New Zealand: Birch Hill (Lake Pukaki basin) (Putnam et al., 2010a); Irishman Stream (Lake Ōhau basin) (Kaplan et al., 2010) and Whale Stream (Lake Pukaki basin)  (based on two different moraine ridgestwo ages from per moraine). (e) Normal kernel density plots from South America: Torres del Paine-II (García et al., 2012); Río Tranquilo-4 (Sagredo et al., 2018); Lacteo lateral (Mendelová et al., 2020). (g) Normal kernel density plots from the upper right tributary of the Ahuriri River valley (current study). [Color figure can be viewed at wileyonlinelibrary.com] variables because of their control on the elevation of the snow/ rain threshold, and also the melt process via the exchange of turbulent heat . In addition, glacier modelling studies show that in effect this means that large precipitation changes (much larger than what we see in present-day interannual variability) are needed to offset the effects of relatively small temperature changes on New Zealand glaciers (Anderson and Mackintosh, 2006). For example, Doughty et al. (2013) showed that a c. 3.5× increase in precipitation (relative to modern levels) is required for precipitation alone to drive glacier advance during the ACR. This study and that of Eaves et al. (2017) also show that more reasonable precipitation changes of ±20% only alter Lateglacial temperature estimates by <0.5°C. Thus, while it is likely that precipitation may have varied during the Lateglacial (e.g. Whittaker et al., 2011), the absence of quantitative estimates, coupled with the relatively low sensitivity of New Zealand glaciers to this variable, means that we choose to present our ELA-based climate estimates purely as a temperature signal. The studies above and the conclusions of Rowan et al. (2014) indicate that this assumption imparts uncertainty to our temperature reconstructions on the order of <0.5°C.
Our estimate of ≤700 m snowline depression at 14.5 ± 0.3 ka based on the MELM method equates to air temperatures ≤3.9°C lower than today , assuming no changes in past precipitation. This finding compares to the previously estimated LGM and post-LGM ELA depression from the Ahuriri River valley by ∼880 and ∼770 m at 19.8 ± 0.3 and 16.7 ± 0.3 ka, respectively (Fig. 12) (Tielidze et al., 2022). These previously reconstructed ELA anomalies implied that local air temperature was 5.0 ± 1.0°C lower than present  at 19.8 ± 0.3 ka, while it was 4.4 ± 0.9°C lower at 16.7 ± 0.3 ka, assuming no change in precipitation (Tielidze et al., 2022). Together these results suggest there was at least 0.5°C of warming between 16.7 ± 0.3 and 14.5 ± 0.3 ka (Fig. 12).
Our MELM-and AAR-based estimate of snowline depression by ≤630 m and~360 m along with ≤3.5 and 2.3 ± 0.7°C cooler climate at 13.6 ± 0.3 and 12.6 ± 0.2 ka from the upper right tributary of the Ahuriri River valley is comparable with other Lateglacial studies from New Zealand. For example, AAR analysis of reconstructed glaciers by Porter (1975) show that the Lateglacial ELA in the Tasman River catchment (inside the Birch Hill moraine limit) was 500 ± 50 m lower than modern. Numerical modelling results by Kaplan et al. (2013) suggest that a mean annual temperature of 2.2 ± 0.4°C lower than today is required to simulate a glacier to the Lateglacial (between 15.4 and 12.9 ka) moraine loop at the Whale Stream area (Pukaki basin). Manual glacier reconstruction and glacier modelling suggest that the ELA of the former Otira Glacier near Arthur's Pass was 540-330 m lower between 16 and 14 ka, corresponding to a 3.5-2.2°C colder climate than present-day (Eaves et al., 2017). The coupled energy-balance and ice-flow model also demonstrate an annual temperature~3.2-2.3°C lower at 13.0 ka than present, at nearby Irishman Stream, corresponding 400 m lower ELA than present-day   (Fig. 12). Numerical glacier modelling from Mt. Ruapehu in the central North Island of New Zealand show that the ELA on Mt. Ruapehu was c. 400 m lower than present between 14 and 11 ka, which equates to a temperature anomaly of 3-2°C lower than present (Eaves et al., 2019). Our snowline-based temperature increase of~1.2°C between 13.6 ± 0.3 and 12.6 ± 0.2 ka from the upper right tributary of the Ahuriri River valley is also comparable with a temperature increase of 0.25-1.0°C from the Irishman Basin between 13 and 12 ka (Kaplan et al., 2010).
Multiple atmospheric and oceanic climate proxy records from across Southern high latitudes reveal that the cooling of the Lateglacial period started at around 14.6 ka and reached its cooling peak at the end of the ACR after~13 ka (Fig. 11a-d). This is in contrast to our results which suggest lower temperatures at 14.5 ± 0.3 and 13.6 ± 0.3 ka than at 12.6 ± 0.2 ka. Our suggestion of lower temperatures early in the Lateglacial is more similar to the radiocarbon-dated pollen and chironomid record from Boundary Stream Tarn at Lake Pukaki (Fig. 13a). These data suggest that early cooling of the Lateglacial occurred between 14.2 and 14.0 ka when temperatures decreased to 1.5-2.0°C below modern summer  air temperature (Vandergoes, et al., 2008). The onset of this cooling period at 14.2 ka matches within uncertainty with our oldest moraine ages of 14.5 ± 0.3 ka. A short interruption to this cooling occurred at~13.9 ka, when mean summer temperature increased to around 2.5°C, or slightly above present-day level. Over the next 300-500 years (i.e. between 13.6-13.4 ka), average mean summer temperatures dropped to around 2.0-2.9°C below the present-day and in the most extreme case as much as 3.9°C lower (Vandergoes, et al., 2008) (Fig. 13a). Again, this fits within uncertainty of our moraine ages of 13.6 ± 0.3 ka. Isotope ratio profiles from Hollywood Cave stalagmite (HW3), South Island, New Zealand, also suggest that the coldest Lateglacial climate occurred between 14.2 and 13.9 ka (Whittaker et al., 2011) ( Fig. 13b) coinciding within uncertainty of our 14.5 ± 0.3 ka ages. Furthermore, five different ice core records from Antarctica (Law Dome, Byrd, Siple Dome, Dronning Maud Land and Talos Dome) show a strong cooling trend from the onset of the ACR at c.~14.6 ± 0.2 to~14 ± 0.2 ka, then a period of variability without significant cooling until~13 ka (Pedro et al., 2011;Stenni et al., 2011) (Fig. 13c-g), again consistent with our lowest temperatures from outer moraine ages of 14.5 ± 0.3 ka.

Possible causes of early cooling
The ACR coincides with the Bølling-Allerød warm stage in the Northern Hemisphere, which is thought to be an example of the inter-hemispheric coupling of abrupt climate change (e.g. Stocker and Johnsen, 2003;Barker et al., 2009;Newnham et al., 2012;Pedro et al., 2016). Stocker and Johnsen (2003) suggested that the gradual cooling in Antarctica during Greenland interstadials such as the Bølling-Allerød results from enhanced cross-equatorial northward ocean heat transport by the Atlantic Meridional Overturning Circulation (AMOC). In their view, increased northward heat transport during Greenland interstadials abruptly cools the South Atlantic, and more gradually cools the Southern Ocean, due to its larger thermal inertia. More recent work attributes the Southern Ocean and Antarctic cooling to gradual mixing of the South Atlantic cold signal across the Antarctic Circumpolar Current (ACC) along with sea ice feedbacks (Schmittner et al., 2003;Pedro et al., 2018). Rapid ocean propagation of temperature anomalies is also possible (north of the ACC) by Kelvin and Rossby wave-driven adjustment of isopycnals and thermocline depths (Pedro et al., 2018). Additionally, there is growing evidence for a fast atmospheric teleconnection operating at the onset of the ACR, and indeed at the onset of all millennial-scale Antarctic cold events during the last glacial (Markle et al., 2017;Buizert et al., 2018). Specifically, ice-core deuterium excess records indicate abrupt northward shifts of the moisture source regions at Antarctic ice core sites within decades of the onsets of Greenland interstadials throughout the last glacial period, which is consistent with northward-shifted westerlies. Models suggest that atmospheric adjustments during these millennial-scale events are a thermodynamic response to a change in the energy balance between the hemispheres during abrupt climate change. At the onset of Greenland interstadials, the abrupt northern warming shifts the thermal equator northward (Stocker and Johnsen, 2003). The atmosphere responds by fluxing more heat and momentum southward, toward the cooler hemisphere, in the crossequatorial Hadley circulation. The increased momentum flux into the Southern Hemisphere draws the westerlies north-  Figure 13. (a) Boundary Stream Tarn (Lake Pukaki area) chironomid temperature records (Vandergoes, et al., 2008). (b) Isotope ratio profiles (five-point means) from stalagmite HW3 from Hollywood Cave, South Island, New Zealand (Whittaker et al., 2011). (c-f) Deglacial δ 18 O ice core records from Law Dome, Byrd, Siple Dome and Dronning Maud Land (DML), Antarctica (Pedro et al., 2011). (g) Deglacial δ 18 O ice core record from Talos Dome, Antarctica (Stenni et al., 2011). [Color figure can be viewed at wileyonlinelibrary.com] ward (Ceppi et al., 2013;Pedro et al., 2016). The thermal imbalance between the hemispheres is largest within the first centuries of the Greenland interstadials and thereafter relaxes as the interstadial progresses (Clement and Peterson, 2008). On this basis, we propose that the early cooling and maximum Lateglacial ice extent in New Zealand in the first centuries of the ACR reflects the northward shift of the westerlies, causing lower temperatures and/or relatively higher precipitation. Subsequent gradual southward relaxation of the winds, causing warming and/or precipitation reduction, resulted in gradual glacier retreat as seen across the Lateglacial moraine records in the Southern Alps (this study; Kaplan et al., 2013). These wind changes may have been accompanied by sea surface temperature changes north of the ACC induced by ocean propagation of abrupt climate change. Discriminating between the relative roles of the atmosphere and ocean is difficult in this highly coupled system. However, we argue that the atmospheric component provides a more physically consistent explanation for the pattern of glacier length changes emerging in southern mid-latitude terrestrial records.

Conclusions
In this paper, we aimed to improve our understanding of the glacier behaviour and climate events during the ACR in the Southern Alps of New Zealand using cosmogenic exposure dating and snowline reconstructions. Based on several moraine features and 28 10 Be ages, we reconstructed the geometries of former glaciers in two different tributary river basins of the Ahuriri River valley. Using the AAR and MELM methods, we also estimated associated snowline elevations. Furthermore, based on the reconstructed glacier snowline we estimated palaeo temperatures during the ACR. Our results provide insight into the structure of the ACR event in New Zealand. The 14.5 ± 0.3 ka outer moraine feature from the upper right tributary of the Ahuriri River valley confirms the largest glacier advance in New Zealand occurred in the first centuries of this millennial-scale event. By 12.6 ka, the local snowline had risen~340 m (from~1300 to~1640 m a.s.l.), which corresponds to~1.6°C of warming during the ACR. We argue that this pattern of early glacier advance followed by gradual retreat could be explained by an abrupt northward shift of the southern westerly winds in the first centuries of the ACR, followed by gradual southward relaxation during the remainder of this southern stadial event.

Data availability statement
The data that supports the findings of this study are available in the supplementary material of this article.

Supporting information
Additional supporting information can be found in the online version of this article. Table S1. The exact data block that was entered into the calculator for outer moraine ridge (44°7′57″ S 169°38′34″ E). Table S2. The exact data block that was entered into the calculator for inner moraine ridge (44°7′55″ S 169°38′30″ E). Table S3. The exact data block that was entered into the calculator for terminal moraine ridge (44°7′54″S 169°38′42″E). Table S4. The exact data block that was entered into the calculator for the Riegel-bedrock near the inner terminal moraine ridge (44°7′56″S 169°38′44″E). Table S5. The exact data block that was entered into the calculator for Canyon Creek bedrock (44°11'22″S 169°35'28″E).