Lake sediments uncover glacier change on the Kerguelen Islands during the past 1500 years

1 Introduction

Anthropogenic warming and shifts in atmospheric circulation are driving the retreat of sub-Antarctic glaciers (Gordon et al., 2008; Favier et al., 2016; Nel et al., 2023). As sensitive climate indicators (Haeberli et al., 2007), glaciers respond rapidly to environmental change, making palaeo-glacier reconstructions a valuable tool for placing current climate trends in a longer-term context. Such reconstructions can assess whether recent retreat exceeds natural variability, identify underlying climate drivers, and ultimately improve predictions of future glacier-climate behaviour (Solomina et al., 2015). Here, we aim at reconstructing past glacier variability on the sub-Antarctic Kerguelen Islands (49°S), located in the Indian Southern Ocean (Figure 1), based on sediments from two glacier-fed lakes (Figure 2b). Such lakes can preserve continous archives of glacier activity in the catchment (Karlén, 1976), created by the delivery and accumulation of fine-grained sediments produced by glacial erosion (Leemann and Niessen, 1994; Hallet et al., 1996). Because the amount of sediment transported each melt season is proportional to glacier extent in temperate systems (Dahl et al., 2003) like those found on the Kerguelen Islands, the accumulation of glacial sediments in glacier-fed lakes can serve as a proxy for past (relative) glacier fluctuations (Karlén, 1981; Bakke et al., 2009). Using this approach, we here aim at reconstructing glacier activity on the northeast side of the Cook Ice Cap (403 km²) (Cogley et al., 2014), the largest ice mass on the Kerguelen Islands.

Southern Ocean climate is dominated by the Southern Hemisphere westerlies (SHW), a zonal wind belt that moves around Antarctica (Figure 1). Changes in the intensity and latitudinal position of the SHW is modulated by the Southern Annular Mode (SAM), the dominant mode of atmospheric climate variations in the extratropical Southern Hemisphere (SH) (Gong and Wang, 1998; Fogt and Marshall, 2020). The SAM operates in two distinct phases: (i) a positive phase in which the SHW move poleward and gain strength, and (ii) a negative phase in which the SHW expand equatorward and weaken (Fogt and Marshall, 2020). Alternations between the SAM modes significantly impact temperature, precipitation patterns and currents across the Southern Ocean (Biastoch et al., 2009; Le Quéré et al., 2007; Gillett et al., 2006).

Over the past 70 years, the SAM has remained positive, more so during the austral summer (Fogt and Marshall, 2020). This trend could have resulted in accelerated Antarctic Peninsula warming (Vaughan et al., 2003) and extreme weather events in the mid-latitude SH such as, droughts (Cai et al., 2011), wildfires (Holz et al., 2017) and extreme rainfall (Huang et al., 2021). Moreover, the increased wind stress from the intensified SHW have strengthened the Antarctic Circumpolar Current (ACC), hence weakening the ability of the Southern Ocean to sequester atmospheric CO² (Lovenduski et al., 2008; Lenton and Matear, 2007). Climate projections indicate a continuation and possible further escalation of this trend, although estimates differ (Zheng et al., 2013; Goyal et al., 2021; Deng et al., 2022; Le Quéré et al., 2007). Given the uncertain future of the SAM and its central role in regulating Southern Ocean (and by extension global) climate, understanding its past (and natural) variability is essential, especially considering recent SAM trends are amplified by anthropogenic drivers such as greenhouse gas-induced warming and the loss of stratospheric ozone (Thompson and Solomon, 2002).

Terrestrial proxy records from regions influenced by SAM dynamics, such as the Indian Southern Ocean, can help extend our knowledge baseline of this system (Fogt and Marshall, 2020). Based on observations, there is indeed a significant relationship between the climate around the Kerguelen Islands and the SAM (Gillett et al., 2006). This is supported by a modelling study, which indicate 77% of the observed retreat at the Cook Ice Cap is due to reduced moisture availability from the persistent positive SAM phase and southward shifted SHW (Favier et al., 2016). We may thus assume that during negative SAM phases the opposite occurs, as the SHW expands equatorward, bringing colder and wetter (snowier) conditions that favour glacier advances. In line with this assumption, previous glacier reconstructions from the archipelago suggest SAM-like variability and sea surface temperature (SST) changes are the primary drivers of glacier fluctuations during the Holocene (Charton et al., 2022; Tuestad et al., 2024; Chassiot et al., 2024; Sleire et al., 2025; Jomelli et al., 2017). The influence of SST is owed to the extremely maritime setting (see 2. Setting). It has hence been shown that over the last 8000 years, mean annual air temperature (MAAT) and SST evolved in a roughly similar way (Bellet et al., 2026). The Kerguelen Islands currently sits in the direct path of the ACC (Figure 1), between the Sub-Antarctic Front (SAF) and Polar Front (PF) (Park et al., 2019). The PF constitute a seasonally migrating water mass boundary (Pauthenet et al., 2018), that influences the ocean temperature around the archipelago. Similar PF migrations are also suggested to operate on Glacial-Interglacial (Civel-Mazens et al., 2021a) as well as Holocene timescales (Bellet et al., 2026). Building on previous Holocene glacier reconstructions (see above), here we date the maximum extent of the Agassiz Glacier (17 km²) and reconstruct continous fluctuations of the Chamonix Glacier (12 km²) (Figure 2b) (Deline et al., 2024), by analyzing sediments from two lakes: Lake Athéna II and Lake Héra (Figure 3), using a multiproxy approach that includes X-Ray Fluorescence (XRF) geochemistry, magnetic susceptibility (MS), dry bulk density (DBD), loss-on-ignition (LOI), and grain-size analysis, combined with radiocarbon (¹⁴C) dating.

2 Setting

The Cook Ice Cap, which feeds the investigated glacier-fed lake basins, is situated on the SHW-exposed western side of the 6700 km² and 150 km wide Grande Terre, the main island of the archipelago. The bedrock in our study area is comprised mainly of basaltic lava flow rocks (Giret et al., 2003) with some light-coloured intrusions (Mathieu et al., 2011). On the drier eastern Kerguelen Islands, at the Courbet Peninsula (Figure 2a), lies the French research base Port-Aux-Français (PAF). PAF weather data span year 1951 Common Era (CE) to the present and reveal an oceanic climate with a MAAT of 4.5 °C and 800 mm of precipitation (Berthier et al., 2009). The Chamonix Glacier flows north-eastward from the Cook Ice Cap. Based on estimates from the adjacent Ampère Glacier located on the southeastern side of the ice cap (Berthier et al., 2009), as well as other similarly-sized maritime glaciers in the mid-latitude SH (Raper and Braithwaite, 2009), the Chamonix Glacier likely has a 25–30 years response time to climatic forcing.

The glacier foreland is characterised by a meandering proglacial river cutting into a massive, dark, grey-coloured basalt outwash plain, dotted with numerous pronounced till ridges that constitute extensive sequences of latero-frontal moraines that have been mapped in detail (Deline et al., 2024). A total of eight sequences of moraines are located between Lake Chamonix and Lake Aphrodite (Figure 2b). Based on cosmogenic exposure ages (³⁶Cl) the outermost moraine ridges ahead of the Chamonix Glacier were deposited 1700 ± 800 years ago, while the third sequence is suggested to have been deposited 320 ± 50 years ago (Tuestad, 2024). During the 2019 PALAS II campaign, long (piston) and shorter (gravity) sediment cores were extracted from deep and flat (undisturbed) sub-basins in echosounder-surveyed (see 3. Methods) glacier-fed lakes (Figures 3, 4). To ensure the longest possible sediment history of past Chamonix Glacier fluctuations, emphasis was placed on cores retrieved from Lake Athéna II (0.35 km²) at 100 m a.s.l. and Lake Héra (1.25 km²) at 80 m a.s.l. (Figure 3).

Lake Athéna II is connected to the more proximal main Lake Athéna through a shallow, 2-m wide sill. From here, a proglacial river carries the remaining suspended sediment load from the Chamonix Glacier to the distalmost lake; Lake Héra. Reconstructions of the glacier margin of the neighbouring Agassiz Glacier (Deline et al., 2024; Tuestad, 2024) (Figure 2b) suggest that during the deposition of the outermost moraines the ice front extended past the current outlet of Lake Héra, likely allowing for the direct input of meltwater and glacier-derived material. As such, Lake Héra has likely functioned as a glacial threshold lake, only receiving Agassiz Glacier sediments when the ice front exceeds a topographic boundary (Larocca and Axford, 2022). While these lakes cannot resolve continous changes, they can be used to accurately date past glacier extents (Briner et al., 2010; Håkansson et al., 2014). Cosmogenic (³⁶Cl) exposure ages from the outermost moraine ridges loosely constrains the age of deposition to between 2900 and 5000 years ago (Tuestad, 2024). Subsequent glacier advances that led to the deposition of the other moraines were not significant enough for the Agassiz Glacier to deliver sediments to the lake.

3 Methods

3.1 Coring and Quaternary geological mapping

Well-suited sub-basins for coring were identified by first mapping the lakebed bathymetry and the thickness and structure of the soft-sediment infill in the lakes. For this, we relied on a Lowrance Elite TI² sonar for the bathymetry and a EDGETECH 3100 Compressed High-intensity Radar Pulse (CHIRP) for the soft-sediment infill. The CHIRP transmits a sweeping range of frequencies (2–16 kHz). The CHIRP data were later processed and visualised using SeiSee (v. 2.2) (Figure 4), where a band-pass filter was applied and the colour contrast of the reflectors were enhanced. Two-way travel time was converted to depth using an average velocity of 1500 m/s. The sonar depth measurements used to reconstruct the bathymetry were automatically assigned coordinates during sampling. In ArcGIS Pro (v. 2.9.5), digitized lake outlines were converted to point data and combined with the sonar depth measurements, yielding approximately 150.000 data points for Lake Athéna and Lake Athéna II combined and approximately 24.000 data points for Lake Héra. Interpolating between these data points using the Topo to Raster tool generated the lakebed bathymetries shown in Figure 3. Following the lake survey, we extracted pairs comprised of a long piston core and a shorter gravity core to preserve the water-sediment interface. For this, we used two different coring devices: a customized piston corer; a Nesje corer (Nesje et al., 1987) and a UWITEC gravity corer (see Table 1). The piston corer consists of a 110 mm PVC liner equipped with a greased piston and core catcher. The liner was lowered down to the lakebed and pushed into the sediments by manually raising and dropping a weight using ropes. The UWITEC gravity corer consists of a 64 mm liner without a core catcher to more easily preserve the uppermost (and most easily disturbed) layers of sediment. To link lake sediments to different catchment processes, the geomorphology around the study lakes was mapped. The field-mapping results were later supplemented with drone images (n = 155) of the study lakes and associated landforms and deposits (Supplementary Figure 2) acquired using a DJI Phantom drone. In addition, a high-resolution (0.5 m) satellite image (DS_PHR1A_201901160513269_FR1_PX_E069S50_0121_01703/696/690) was used to identify larger-scale geomorphic features across the study area.

TABLE 1. Main characteristics of the collected sediment cores analysed in this study.

3.2 Laboratory analysis

We focused our analysis on the core pairs HEP-19–02/HEG-19–02 extracted from Lake Héra and AT2-19–04/AT2-19–05 from the upstream Lake Athéna II (Figure 3) based on core lengths (see Table 1) and the potential (or lack thereof) for disturbance and/or dilution of the glacial sediment signal (see 4.1 Mapping and lake survey). Furthermore, selected analysis was performed on a second Lake Héra core pair; HEP-19–03/HEG-19–03, to aid in constructing a master record from that lake. The cores were processed at the EarthLab sediment laboratory in Bergen. First, the cores were split lengthwise and cleaned. Afterwards, working and reference halves were designated and stored in a temperature-controlled (∼4 °C) setting. Using the ITRAX XRF scanner from Cox Analytics (Croudace et al., 2006) the working halves of the cores were photographed and the elemental geochemistry was mapped at 0.02–0.05 cm resolution (AT2-19–04 n = 10,116, AT2-19–05 n = 1750, HEP-19–03 n = 3580, HEG-19–03 n = 1981, HEP-19–02 n = 9627, HEG-19–02 n = 4274). The scanner was fitted with a Molybdenum (Mo) tube, well-suited for detecting heavier elements, e.g., Calcium (Ca) and Titanium (Ti) (Croudace, 2015) that can dominate glaciolacustrine sediments (Bakke et al., 2009), and utilised power settings of 35 kV and 30 mA with a 10 s dwell-time.

To check for missing surface sediment Ti counts per second (cps) from the piston- and gravity core pairs were aligned on a common depth scale (Supplementary Figures 1d,e) in the version 9.10.0 R2021A of the MatLab software (MathWorks, 2021) following the approach by Trauth (2020) that involves linearly interpolating between tie-points. To normalise the relative XRF values and account for non-linear matrix effects, e.g., water content, elemental counts of Silicon (Si), Calcium (Ca), Titanium (Ti), Manganese (Mn) and Iron (Fe) were centred log-ratio (clr) transformed (Bertrand et al., 2024). Next, a Multisensor Core Logger was used to measure surface MS at 0.5 cm (AT2-19–04 n = 1011, HEP-19–03 n = 718, HEG-19–03 n = 198, HEP-19–02 n = 972, HEG-19–02 n = 172). Subsequently, the clr-transformed elements were used in a Principal Component Analysis (PCA) using XLSTAT (Addinsoft, 2019) to compute the main trends (PCs) in the XRF datasets with MS plotted on the ordination diagrams as a supplementary variable. Computed tomography (CT) was used to assess relative density changes in the transition between units 2 and 1 in HEP-19–02 composite (see 4.2 Stratigraphy). Using the ProCon CT-ALPHA scanner helical scans were generated. Tube voltage, current and exposure time were set to 130 kV, 850 mA and 334 m, producing 1600 images. Beam-hardening effects were reduced by filtering with a 1 mm Al-disc (Brooks and Chiro, 1976). Post-processing was performed using version 9.1.1 of the ThermoFisher Amira-Avizo software to generate a 2-D CT image (XY; ortho slice) of the scanned interval. Furthermore, the Volume Rendering module was used to visualise larger clasts in 3-D.

To identify downcore trends in physical sediment properties, we sampled 1 cm³ cylinders of sediment every centimetre throughout HEP-19–02 (n = 225). Dry bulk density (DBD) and water content were measured after drying at 50 °C for 48 h, and organic content was approximated by Loss-On-Ignition at 550 °C (LOI⁵⁵⁰) for 1 h (Dean, 1974). To measure downcore trends in grain-size distributions, every 20th sample was analysed following incineration (Vasskog et al., 2016; van der Bilt et al., 2021). For this, we relied on the Malvern Panalytical Mastersizer 3000 connected to a Hydro LV wet-dispersant unit running at 2500 rotations per minute. Sodium hexametaphosphate and an ultrasound treatment (60 s) were used to disperse individual particles before analysis. Each sample was analysed in quintuplets before the average was calculated and used as the final measurement. Grain-size statistics were calculated using the GRADISTAT spreadsheet (Blott and Pye, 2001).

3.3 AMS radiocarbon dating and short-lived radionuclide activities

Plant remains were wet sieved at 250 µm to remove adhering sediments. The extracted samples (AT2-19–04 n = 10, HEP-19–03 n = 6 and HEP-19–02 n = 8) were dried at 50 °C for 48 h, weighed and placed in sterilised vials to be submitted for Accelerator Mass Spectrometry (AMS) radiocarbon dating (¹⁴C) at the Poznan Radiocarbon Laboratory (Poz) and Lund Radiocarbon Dating Laboratory (LuS). To constrain surface sediment age, bulk material (>2 g dry weight) was extracted every 0.5 cm from HEG-19–02 (n = 15) and short-lived radionuclide (¹³⁷Cs, ²⁴¹Am, ²¹⁰Pb, ²²⁶Ra) activities were established using the EDYTEM SAGe well-type detector (Reyss et al., 1995) installed at the Laboratoire Souterrain de Modane (LSM) underground laboratory. The total ²²⁶Ra activity was subtracted from the total ²¹⁰Pb activity to calculate excess ²¹⁰Pb (²¹⁰Pbexc) (22.3 years half-life). Regrettably, due to low activity of short-lived radionuclides in the uppermost samples (Supplementary Figure 1b), it was not possible to model an age-depth relationship. Therefore, the adjacent gravity core HEG-19–03 was used to capture the most recently deposited sediments in Lake Héra (see 4.2 Stratigraphy).

4 Results

4.1 Mapping and lake survey

As highlighted above (see 2. Setting), a previous study has mapped the moraines in the Chamonix Glacier foreland, as well as the neighbouring Agassiz Glacier (Deline et al., 2024). Here, we focus on landforms and deposits around our study lakes to get an overview of potential processes within their respective catchments, a prerequisite for an accurate interpretation of the sediment in the lakes (Bakke et al., 2010). Towards the northeast, Lake Athéna II is surrounded by gentle slopes (∼1°–2°) with little visible vegetation. To the west are steeper areas (∼20°) with multiple alluvial fans of varying size, between 700 and 3800 m², likely formed by streams reinvigorated by seasonal snow melting. To the east is a 100 m long beach with an 11 m shoreline sequence. Cutting into the beach is an active stream that originates from the nearby mountain. Based on the lake survey, Lake Athéna II’s maximum depth is 17 m, with the coring location for AT2-19–04/AT2-19–05 at a depth of 14 m, respectively. The core pair was retrieved around 120 m from the beach. Seismic profile (SP) 1 depicted in Figure 4a shows the seismostratigraphy of the coring location, which we interpret as two principal acoustic facies. Facies 1 consists of laterally continous, high-amplitude, parallel to subparallel reflections forming a draping unit that follows the underlying basin morphology. The setting of the lake and the draping indicate deposition of suspension in a glaciolacustrine environment following basin deglaciation. Facies 2 starts at around 25 m depth. It’s a mostly acoustically transparent interval with a gently undulating upper surface. Internal reflectors are largely absent, likely owing to poor acoustic penetration into coarse, consolidated material. The position of this facies beneath the layered sediment succession of Facies 1, together with its seismic character, suggests it represents glacial till and/or bedrock.

Moving downstream, Lake Héra sits between steeper (∼23° to the northeast and ∼13° to the southwest) slopes mostly covered with till. The northern tip of the lake is characterised by a 0.11 km² fluvial delta with well-developed, braided rivers. Next to the delta is colluvium with boulders deposited along the edge of the lake. A similar beach deposit (∼85 m long) to that observed at Lake Athéna II is found towards the east, where the proglacial river enters the lake (Supplementary Figure 2b). Distinct moraine ridges from the Agassiz Glacier are present at the lakes’ outlet, where mulitple relict meltwater channels are carved into glaciofluvial outwash sediments. Combined, these features suggest that meltwater from the Agassiz Glacier spilt over into this lake in the past (Figure 3) (Deline et al., 2024). The lake bathymetry is complex with multiple sub-basins. Based on the bathymetric map (Figure 3), the main basin reaches a maximum depth of ∼62 m, respectively. The core pairs HEP-19–02/HEG-19–02 and HEP-19–03/HEG-19–03 were extracted from depths of 7 m and 31 m, respectively. The main basin was avoided during coring as SP 2 reveals locally disturbed sediments characterised by discontinuous and disrupted reflectors. The disturbance of the main basin could be related to subaqueous mass wasting as indicated in Figure 4b. Despite these local disturbances, SP 2 and SP 3 indicate the same two principal acoustic facies as observed in Lake Athéna II, although greater lateral variations in facies thickness and reflection amplitude. Around the HEP-19-02/HEG-19-02 coring site (SP 3), Facies 1 forms a ∼6 m thick draping unit composed of high-amplitude, parallel reflectors that transitions into an approximately 8 m thick acoustically semi-transparent interval with fewer discernible internal reflections. At the HEP-19-03/HEG-19-03 coring site (SP 2), Facies 1 is approximately 8 m thick and is characterised by continuous, high-amplitude, horizontal reflectors that extend laterally across the coring site (Figure 4b). The recognition of the same acoustic facies in the downstream Lake Héra supports our interpretation that Facies 1 is a glaciolacustrine unit while Facies 2 is glacial till and/or bedrock, although limited penetration prevents a clear distinction between these two.

4.2 Stratigraphy

The primary piston cores AT2-19–04 and HEP-19–02 were aligned with gravity cores (AT2-19–05 and HEG-19–03) (Supplementary Figures 1c,d) to ensure preservation of the surface sediments by linearly interpolating between tie-points (n = 5, 8) (Trauth, 2020) (see 3. Methods). Because of missing surface sediments in HEG-19–02 (see 4.3 Geochronology), the gravity core HEG-19–03 was used to capture the most recently deposited sediments in Lake Héra, revealing a 51 cm gap at the top of the HEP-19–02 piston core (Supplementary Figure 1e). At the same time, close agreement between AT2-19–04 and associated gravity core AT2-19-05 suggests negligible loss of surface sediments (Supplementary Figure 1c) in the AT2-19–04 piston core. However, comparing the core top of AT2-19–05 with HEG-19–03 suggests significant sediments are missing from the Lake Athéna II cores. Sediments in core AT2-19–04 and the HEP-19–02 composite record were divided into units and subunits based on XRF and MS signatures. Units and subunits that are correlated between the records, e.g., Subunit 1b (Figure 5), represent the same type of sediments presumed to be deposited at the same time, and under the same environmental conditions.

Basal Unit 2 (240 cm), found only in the HEP-19–02 composite record (Figure 5b), is a dark, grey-coloured sediment intercalated with thick light, grey-coloured lenses. Based on analysed grain-size distributions (n = 15), the sediment falls in the (very) coarse silt fraction (Figure 5c), and consistently high DBD values (ca. 1.95 g/cm³) reveal the sequence is dense. Moreover, high values in clr-transformed Titanium-clr (Ti-clr) reflect a clastic character (Bakke et al., 2009; Haug et al., 2003). This is demonstrated statistically in the ordination diagram (Figure 7b) that show Unit 2 samples cluster near Ti-clr, as well as Ca-clr, Si-clr and supplementary variable MS. MS is a proxy for magnetic mineral abundance in sediments, a parameter that is closely linked to the overall clastic content (Thompson et al., 1975; Støren et al., 2008). This clastic characterization is further evidenced by reduced organic content in Unit 2 detected by (i) low LOI⁵⁵⁰ percentages and (ii) positively correlated Incoherent-to-Coherent (Inc./Coh.) scatter ratios, a readily applied XRF ratio for estimating the organic component of sediments (Liu et al., 2013). The CT image in Figure 6a reveals an abrupt shift from dense, uniform sediments (Unit 2) to mm-scale laminated silty bands (Unit 1, Subunit 1e). Furthermore, the CT data shows semi-rounded pebbles nested in the bottommost laminations of Unit 1.

Fine-textured beige-brown-coloured Unit 1 sediments exhibit greater variability in the analysed proxies. Based on this, we further subdivided Unit 1 in both records into subunits 1a – 1e, with the latter only found in HEP-19–02 composite (Figure 5b). Subunits 1e, 1c (152 cm in AT2-19–04 and 51 cm in HEP-19–02 composite) and 1a (257 cm in AT2-19–04 and 128 cm in HEP-19–02 composite) display high values in Ti-clr and MS. Moreover, stable and low LOI⁵⁵⁰% and Inc./Coh. scatter ratios (Figure 5) both indicate a low concentration of organic matter in the sediments. Dry bulk densities in HEP-19–02 composite samples reveal these units are, on average, denser than the more recent and more organic-rich ones (Figure 5b), in coherence with the predominance of very fine silts in the grain-size spectra.

The opposite trend is observed in Subunit 1d (47 cm in AT2-19–04 and 28 cm in HEP-19–02 composite), which is notably coarser in texture (Figure 5c). As a result, the sediments are less densely packed (x ¯ DBD = 0.76 g/cm³). The organic content is high (x¯ LOI⁵⁵⁰ = 3.4%), and we note lows in clastic indicator MS (x¯ Si 10⁻⁵ = 1265 in AT2-19–04, x¯ Si 10⁻⁵ = 647 in HEP-19–02 composite). Furthermore, Subunit 1d (in AT2-19–04 only) (Figure 5a) includes a dark-coloured layer with tephra shards. Subunit 1b (48 cm in AT2-19–04 and 14 cm in HEP-19–02 composite) presents the same organic fingerprint as detected in Subunit 1d, culminating with a distinct two-peaked spike in Inc./Coh. scatter ratios present in both records. Furthermore, visual logging of the HEP-19–02 composite record reveals plant residue in the core surface. Simultaneously, here we identify near-identical grain-size distributions to clastic subunits 1a, 1c and 1e. In the ordination diagrams (Figure 7) samples from these subunits (1d and 1b) plot closer to the elements Fe-clr and Mn-clr, and away from Ti-clr, Ca-clr, Si-clr and MS, respectively.

4.3 Geochronology

Core chronologies are based on radiocarbon (¹⁴C) ages derived from plant remains dated at the Poznan (Poz) and Lund (LuS) radiocarbon laboratories. ¹⁴C ages were calibrated using the IntCal20 SH curve (Hogg et al., 2020). We combined terrestrial and aquatic plants in our samples owing to (i) a scarcity of macroscale organic matter fragments and (ii) the lack of carbonates in the catchment area, which leads us to assume the absence of reservoir effects (Zhou et al., 2022). Age-depth relationships were modelled using a Bayesian statistical framework in the R-based (Team, 2021) package Bacon (Blaauw and Christen, 2011). As highlighted above (see 4.2 Stratigraphy), it is likely sediments are missing from the top of AT2-19–04. We therefore conservatively only model an age-depth relationship until the uppermost ¹⁴C age (Poz-147728) (see Table 2). For the AT2-19–04 age-depth model, we identify multiple outlier ages. Anomalously old ¹⁴C ages that falls outside the 2σ (95%) confidence interval are suggested to be caused by re-mobilisation of previously deposited sediments with an abundance of old organic matter fragments (Björck and Wohlfarth, 2001). Regarding the tephra layer in Subunit 1e (see 4.2 Stratigraphy), we regard this section as instantaneously deposited on timescales at which glacier-climate fluctuations take place (Støren et al., 2008) and remove it from the age-depth model using the slump function.

TABLE 2. Overview of radiocarbon (¹⁴C) samples. Note that the depth for samples extracted from HEP-19–02 and HEP-19–03 are modified for the master depth scale (HEP-19–02 composite). Outliers are marked with *.

For the HEP-19–02 composite record, ¹⁴C ages from piston cores HEP-19–02 and HEP-19–03 were combined, and their depths modified to match the master (HEP-19–02 composite) depth scale that includes surface sediments (see above). The near-zero short-lived radionuclide activity (¹³⁷Cs, ²⁴¹Am, ²¹⁰Pb) in the uppermost sediment slice in HEG-19–02 (Supplementary Figure 1c) indicates missing surface sediments. Minimal ¹³⁷Cs traces (5.8 mBq/g^-1) in the core top may suggest deposition before the delayed fallout following the climax of nuclear weapon detonations observed in SH records around 1964–65 CE (Foucher et al., 2021), warranting the use of the adjacent gravity core HEG-19–03 (Figure 3). To test the validity of combining ¹⁴C ages from two sediment cores, as well as aid the process, we placed the investigated Lake Héra piston cores on a common depth scale (Supplementary Figure 1b) using tie-points (n = 10) and the same approach as for matching piston- and gravity core pairs (see 3. Methods). This enabled statistical evaluation of the correspondence between the records, with a mean correlation coefficient (r) of 0.5 (Supplementary Figure 1b). As a result, the final age-depth model is based on thirteen ¹⁴C ages. Using the boundary function, the model was constrained by the stratigraphic transition between units 2 and 1 (Figure 8), which is characterised by a marked change in grain size and density (see 4.2 Stratigraphy). We interpret this transition as an abrupt change in the lake’s boundary conditions, marked by changes in sedimentation rate, thereby supporting the exclusion of samples Poz-183998 and Poz-183009 (see Table 2). However, it is clear that the HEP-19–02 composite chronology has intervals with a relatively high age uncertainty. For example, at 135 cm depth (Subunit 1b) there is a cluster of ¹⁴C ages with considerable differences in their respective age distribution with the 2σ minimum age of the oldest sample being 1158 CE, while the 2σ maximum age of the youngest sample is 1456 CE, respectively.

5 Discussion

5.1 Late Holocene glacier maximum of the Agassiz Glacier

Basal Unit two covers a short timeframe, likely shorter than the calibrated uncertainties (∼85 years) of our chronology. The sedimentology suggests continuous, but high-energy sedimentation of dense clastic material (see 4.2 Stratigraphy). Glacier-derived sediments are usually denser and less organic than authigenic sedimentation, indeed, in this context, in the vicinity of a formerly large glacier (Deline et al., 2024) higher (lower) DBD (LOI⁵⁵⁰) values can reflect added input from a glacial source (Bakke et al., 2005). Combined with geomorphological field evidence (Figure 2): end moraines and relict meltwater channels adjacent to Lake Héra and the coring location (see 2. Setting), we interpret Unit 2 as a sequence of glacial sediments deposited in a distal setting with a short distance between source and sink. This unit, we propose, was therefore likely deposited by the Agassiz Glacier when its front was at or near its maximum extent, delimited by the outermost moraine ridges next to Lake Héra (Figure 2; Supplementary Figure 2). Furthermore, we interpret the semi-rounded pebbles observed in the CT image as lake-ice rafted debris (drop stones). We theorise that these clasts may have been encapsulated, transported, and deposited by lake ice during annual freeze-thaw cycles (Stachowska et al., 2024), following water-level lowering after the cessation of Agassiz Glacier meltwater input. This event would have exposed coarse, fluvially sculpted streambed material initially deposited near the shore of Lake Héra (Medford et al., 2021). A relatively recent and substantial lowering of the water-level is speculative but supported by field evidence: exposed shorelines (see 4.1 Mapping and lake survey; Supplementary Figure 2) to the east with the uppermost ridge raised ∼10 m above the present-day waterline.

We argue, therefore, that the age of the threshold between units 2 and 1 in the HEP-19–02 composite record, modelled using the Bacon.hist function, to between year 349 and 582 CE (2 σ) (x¯ age 495 CE) (Figure 6b), constrains the timing of the Agassiz Glacier’s last retreat from its maximum extent (Figure 9). This age pre-dates the previous estimate for the deposition of the outermost Agassiz Glacier moraine, which is suggested to have occurred sometime between 2900 and 5000 years ago based on two ³⁶Cl exposure ages (Tuestad, 2024), i.e., a minimum of 1500 years prior to the age derived from the sediment record from Lake Héra. However, we note that there are several (3: 1700 ± 500, 2900 ± 400, 5000 ± 400 years ago) conflicting exposure ages for this landform, possibly because the Agassiz Glacier has advanced close to this extent for multiple times during the mid-to late Holocene. Furthermore, we note that the youngest of the ³⁶Cl-derived ages from Tuestad (2024) (1700 ± 500 years ago) aligns with the age for the upper boundary of Unit two in HEP-19–02 composite and suggested cessation of Agassiz Glacier sediment input. Furthermore, this age matches the estimate for the maximum extent of the neighbouring Chamonix Glacier (Tuestad, 2024).

While not the primary focus of this study, this finding supports the idea put forth by Charton et al. (2022) that glaciers on the Kerguelen Islands have behaved differently than the dominant pattern of SH mountain glaciers during the Holocene. The age of the maximum extents of the Ampère Glacier (1110 CE ± 470 years), Arago Glacier (900 CE ± 330 years) (Charton et al., 2022) and Chamonix Glacier (319 CE ± 900 years) (Tuestad, 2024), are all within the last two millennia. On the other hand, glacier reconstructions from other glaciated regions of the SH such as sub-Antarctic South Georgia (Bakke et al., 2021), New Zealand (Putnam et al., 2012), large parts of South America (Reynhout et al., 2019) and references therein, as well as the Antarctic Peninsula (Kaplan et al., 2020) suggest long-term retreat punctuated by progressively smaller glacier re-advances. Charton et al. (2022) associated the relatively late Holocene glacier maximum in the Kerguelen Islands with an increase in SHW-derived moisture and cooling of the nearby ocean, diagnostic features of negative SAM-like conditions in this region.

The only other SH area where this specific pattern of Holocene glacier evolution is observed is southern Patagonia (Reynhout et al., 2019; Hall et al., 2019; Menounos et al., 2013), which also sits in the direct path of the SHW. The correspondence between the glacier evolution in these areas may be because the climate at both locations are highly influenced by SAM-like variability and associated changes in the SHW. However, while we favour the interpretation mentioned above, we cannot exclude the possibility that differences between individual glaciers stem from physical characteristics (e.g., glacier type, size, and/or hypsometry) and/or local factors such as catchment topography (Roe, 2011). Furthermore, different glaciers in the same region, subjected to similar climatic forcings, may respond differently (Røthe et al., 2018). Dated moraines ahead of the nearby cirque glacier Guynemer (Figure 1b), suggests that this glacier reached its maximum extent during the early Holocene and followed the dominant SH glacier trend (Tuestad et al., 2024), demonstrating that this pattern is not ubiquitous across the entire Kerguelen Islands. Our interpretation for the extents of the Agassiz and Chamonix glaciers during the late 5^th century is shown in Figure 9.

5.2 Interpreting glacial sediment flux variations in the context of glacier dynamics

Unit 1 in both Lake Athéna II and Lake Héra consist of fine-grained, dense clastic sediments (Figure 5) deposited during the last ∼1500 years (Figure 8). The PCA of the XRF datasets (see Methods) identifies a dominant geochemical component (PC1) explaining 51.4% and 46.9% of the total variance (Figure 7). PC1 closely tracks MS changes and is characterised by strong positive loadings of Ti-clr, Ca-clr and Si-clr, indicating that changes in lithogenic sediment input is the primary control on sediment composition (Bertrand et al., 2024). Given the setting, we argue the clastic sediment component is predominantly modulated by changes in the input of glaciogenic sediments delivered by the Chamonix Glacier.

Here, we use Ti-clr, a relative value of the intensity of the XRF Ti signal in the sediments, as an indicator of glacial sediment input because it combines high (0.02 cm) sampling resolution with relative insensitivity to redox processes (Bakke et al., 2009). Furthermore, the same parameter and normalisation approach has successfully reconstructed fluctuations of the neighbouring Ampère Glacier (see 2. Setting) using fjord sediments cores from Table Fjord (Chassiot et al., 2024).

The imprint of hillslope erosion/weathering and reworking of older material is difficult to quantity, however, multiple lines of evidence indicate that Ti-clr primarily reflect glacial flour sediments. Firstly, increased Ti-clr signal intensity corresponds with higher proportions of the very fine silt and clay fractions. The influx of these grain-size fractions may be indicative of upstream glacier activity (Larsen et al., 2020). This is supported by a positive correlation between Ti-clr and DBD (r = 0.5, p = <0.0001), because fine-grained sediments have less pore space and, therefore, usually exhibit higher DBD values (Bakke et al., 2005). Furthermore, Ti-clr exhibit lower values during times of added input of coarser, presumably catchment-derived clastic sediments, suggesting local slope erosion and remobilisation exert only a minor influence on the Ti-clr signals. Together, this supports the notion that Ti-clr primarily tracks glaciogenic sediments carried by the meltwater runoff from the Chamonix Glacier.

Having established that Ti-clr reflects glaciogenic input, the remaining question is whether variations in the volume of glacial flour sediments reflect extent changes or other processes. In many temperate glacier-fed lake systems, larger glaciers produce and transport greater volumes of fine sediment because a larger ice-bed interface promotes bedrock abrasion (Dahl et al., 2003). However, this relationship is not necessarily straightforward (Antoniazza and Lane, 2021) and should be evaluated independently for each site. Degradation of cryospheric environments, such as the current retreat documented at the Cook Ice Cap (Verfaillie et al., 2015; Favier et al., 2016), may lead to an increase in meltwater and sediment discharge as previously frozen subglacial tills are thawed and re-mobilised, potentially leading to larger quantities of clastic sediments accumulating in downstream basins despite reduced bedrock erosion rates (Zhang et al., 2022; Delaney and Adhikari, 2020).

Several independent observations suggests that, at our study area, Ti-clr is related to glacier extent. First, Ti-clr shows a gradual decline in the amount of glacial-derived material delivered downstream, mirroring the documented retreat of the Chamonix Glacier (Deline et al., 2024) and Cook Ice Cap as a whole (Verfaillie et al., 2015; Favier et al., 2016). Second, periods of elevated Ti-clr coincide with phases of glacial sediment input in adjacent Lake Aphrodite indicating that the glacier front extended all the way to the threshold south of the lake (Sleire et al., 2025) (Figure 10e). Third, Ti-clr is strongly negatively correlated with LOI⁵⁵⁰ (r = −0.73, p = <0.0001) showing increases in Ti-clr are accompanied by LOI⁵⁵⁰ decreases, consistent with enhanced delivery of glaciogenic sediments by increased bedrock erosion, rather than episodic increases in the meltwater flux. Combined, these lines of evidence suggest that Ti-clr primarily reflects glacier expansion rather than transient increases in meltwater discharge or remobilisation of stored sediments.

We acknowledge several limitations of this approach. We cannot calibrate Ti-clr against instrumental records or dated glacier margins, in large due to the lack of a robust chronology in the most recently deposited sediments, as well as the high chronological uncertainties for the moraine exposure ages (Tuestad, 2024). Likewise, the capacity of the proglacial rivers south of Lake Athéna, and the one connecting Lake Athéna II and Lake Héra, can modify the annual volume of glaciogenic sediments reaching the lakes. Furthermore, glacier erosion rates and consequently sediment availability are also affected by thermal regime changes and what bedrock lithologies are exposed at the sole (Burki et al., 2010). Nevertheless, over the 1500-year period covered by our records, glacier extent has likely been the primary driver of glaciogenic sediment production and delivery over multi-decadal to centennial timescales.

Following from above, we plotted the Ti-clr curves with the mean calibrated ages from the age-depth modelling. However, because parts of the HEP-19–02 composite chronology has a relatively high uncertainty, e.g., Subunit 1b (see 4.3 Geochronology), we aligned it with the AT2-19–04 record using tie-points (n = 15). Duplicate data points were removed to get records with a single data point per year (Trauth, 2020) allowing for a linear regression analysis that shows the aligned Ti-clr timeseries are strongly positively correlated (r = 0.74, p = 0.0001). We preface that the purpose of this is not to over-state the temporal correlation between the two archives: glaciogenic sediments are normally challenging to date accurately due to the clastic character of the sediments and low concentration of plant fragments (van der Bilt et al., 2016), but to more easily evaluate how the glacial sediment signature changes between the two lakes. For example, the more distal record HEP-19–02 composite may be more sensitive to smaller changes in the glacier as less material is transported downstream. AT2-19–04 is favoured for tuning because (i) the chronology is based on ¹⁴C ages derived from a single core and (ii) it has better overlap with reconstructed glacial phases in Lake Aphrodite (Sleire et al., 2025). To avoid interpreting our reconstructions beyond the significant limitations of the chronologies (see above) (Vasskog et al., 2012) as well as smooth out minor Ti-clr changes that may stem from paraglacial reworking (see above), we applied a 5- and 25-year running average smoothing (Figures 10f, g), of which the 25-year smoothed curves forms the basis for the palaeoclimatic discussion (see below).

Our reconstructions of the Chamonix Glacier indicates that the glacier advanced multiple times within the timespan covered by our records with particularly extensive glacier extents between 750 and 1600 CE (Figures 10f, g). This is consistent with local proxy evidence, including a Glycerol Dialkyl Glycerol Tetraether (GDGT)-based MAAT record, which show relatively cold conditions that are interpreted to reflect a predominant negative SAM-like state, between 2000 and 250 cal.yr BP (Figure 10c) (Bellet et al., 2026). Radiolarian- and diatom-based SST records from the Kerguelen Plateau (50°S) and, further south, at the Conrad Rise (54°S), also indicate cool temperatures during this period (Orme et al., 2020; Civel-Mazens et al., 2021b). Cosmogenic exposure dated moraines suggest glaciers across the archipelago advanced during this interval (Jomelli et al., 2017; Verfaillie et al., 2021; Charton et al., 2022; Tuestad et al., 2024; Tuestad, 2024), but the age uncertainties make it difficult to accurately link to changes observed in our records (Figure 10b).

As highlighted previously, the effects of changes in the SAM phase on temperature- and precipitation patterns around the Kerguelen Islands are well-established in observations, models and palaeoclimatic data (see Introduction). Cooling is typically accompanied by atmospheric wetting, while warming coincides with drying, producing strong impacts on local glaciers, as illustrated by the most recent retreat period (Verfaillie et al., 2015; Favier et al., 2016). Given these linkages, glacier advances (retreats), lower (higher) temperatures and increased (decreased) precipitation is expected to co-vary on multi-decadal to centennial timescales during the late Holocene. However, we cannot exclude that the relative significance of different climate variables, and how they correlate varies over longer timescales (Roe, 2011; Trachsel and Nesje, 2015). For example, Bellet et al. (2026) showed a decoupling of air- and ocean temperature variability in the early Holocene. This does not seem to have been the case over the period we consider in the present study (last ∼1500 years) as MAAT, SST and our glacier reconstructions indicate colder conditions along with multiple, extensive glacier advances (see above). Nevertheless, the records presented here suggest there were short-lived periods when the Chamonix Glacier retreated, events that are not obvious from the available temperature reconstructions (see above), likely because of the sampling resolution of these records.

5.3 An apparent 7^th-8^th century CE glacier retreat: reconciling sedimentary evidence with regional palaeoclimate records

As illustrated in Figure 10, the Chamonix Glacier seem to have begun retreating rapidly around 600 CE. The subsequent 150-year interval (Subunit 1d) is marked by a glacier-sediment signal similar to present-day levels, a pattern consistently reproduced across both sediment records. Furthermore, sediments from adjacent Lake Aphrodite suggest a non-glacial system at this stage (Figure 10e) (Sleire et al., 2025). However, the Table Fjord record suggests an advance of the Ampère Glacier during this time (Figure 10d) (Chassiot et al., 2024). Although we cannot exclude that this offset is caused by chronological uncertainties, we consider this unlikely as (i) the Table Fjord records’ chronology is well constrained, and, (ii) while the age-models presented in this study and Sleire et al. (2025) are based on fewer ¹⁴C ages, the retreat at 600 CE is documented at all the investigated lakes, thereby reinforcing the robustness of this finding. It is possible that the behaviour of either glacier were modulated by local conditions, such as whether the fronts were lake-terminating at this time or not, which could dampen/amplify any response to climatic perturbations (Minowa et al., 2015; Harpur et al., 2026). As highlighted by the lack of glacial sediments being deposited in Lake Aphrodite at this stage (Sleire et al., 2025), it is likely the Chamonix Glacier front was indeed lake-terminating within Lake Chamonix. A single data points hints at a slight increase in MAATs around this time (Figure 10c) (Bellet et al., 2026), but there is not sufficient evidence to suggest a pervasive warm period. Moreover, inferred SAM changes in southern Patagonia (51°S), based on non-arboreal pollen, indicates negative SAM-like conditions at this stage (Moreno et al., 2018), at odds with our reconstruction.

The latest SAM index reconstruction, which is modelled using temperature and drought proxy data from sites across the SH, indicates episodic positive austral summer SAM-like phases between 580 and 760 CE (Figure 10a) (King et al., 2023), but the record does not indicate a pervasive positive SAM-like state comparable to present-day observations. Supporting evidence, however, is found elsewhere in areas within the SHW belt. Geochemical- and grain-size data from fjord sediments (Bertrand et al., 2014) as well as changes in pollen ratios in lake sediment records (Moreno et al., 2009; Fletcher et al., 2018), suggest predominantly dry conditions across a broad latitudinal band (46°S - 51°S) of southwestern Patagonia. Latitudinal migrations of the SHW based on lake-level fluctuations (Quade and Kaplan, 2017) and sedimentological and petrophysical analysis (Gilli et al., 2005) at Lake Cardiel in southern South America (49°S), the same latitude as the Kerguelen Islands, suggest the SHW were generally located northwards during the late Holocene, apart from a phase of intensified and presumed southward shifted SHW between 150 and 750 CE. Evidence for a dry period is identified along the Chilen margin (41°S) (Lamy et al., 2001) as well as further northeast towards southern Australia (35°S), where complementary records suggest a dry interval between 500 and 1050 CE linked to positive austral summer SAM-like anomalies (Chamberlayne et al., 2023) and references therein, although a δ ¹⁸O record from ostracod valves from the southeast shows a contradictory pattern of hydroclimate variability (Rahman et al., 2025). At the Falkland Islands (52°S), Turney et al. (2016) suggest intensified SHW, interpreted to reflect drier conditions between 200 and 900 CE. Furthermore, a δ¹⁸O_cell record derived from mosses in a peat bog in southernmost Patagonia (54°S) close to the Drake Passage indicate stronger SHW, expressed as an increase in moisture sourced from the Pacific Ocean between 660 and 810 CE (Xia et al., 2018). At the same latitude and in the same region Waldmann et al. (2010) interpreted a period of increased precipitation reflected in elevated Fe cps in lake sediments as the SHW focused over 54°S, and therefore, potentially, well south of the Kerguelen Islands. As highlighted by Deng et al. (2026) and references therein, opposing SHW behaviour at different localitites may be due to a combined effect of proxy biases, varying age control, local conditions and/or zonally asymmetries in the SHW wind belt. Although inconclusive, it remains a possible explanation that a drier and warmer positive SAM-like climate facilitated the decline in extent of the Chamonix Glacier around 600 CE. Regardless of the forcing mechanism, the glacier-sediment signal produced during this interval implies that the recent retreat documented at the Chamonix Glacier is rare, but not exceptional within the Common Era.

5.4 Short-lived retreats punctuated high glacier activity

Following 750 CE, the Chamonix Glacier re-advanced, with meltwater feeding directly into adjacent Lake Aphrodite (Sleire et al., 2025) (Figure 10e), indicating the glacier front extended beyond Lake Chamonix and, furthermore, was close to the largest extent documented in the moraine record (Figure 1b). The following 850-year interval is characterised by high glacier activity punctuated by several, but short-lived retreat phases, shorter than the chronological uncertainties of our reconstructions. The most pronounced retreat phase that overlaps with multiple of the SAM index reconstructions (Abram et al., 2014; Dätwyler et al., 2018; King et al., 2023) is a two-stage glacier retreat event that occurred between 1150 and 1200 CE that coincided with a sustained positive SAM-like period present in both annual- and austral summer SAM indices (Figure 10a). The Table Fjord record also suggests retreat during this interval, although potential contamination from a layer of tephra makes it difficult to precisely date the retreat phase, as well as the subsequent glacier advance (Chassiot et al., 2024). ¹⁴C dated moss layers ahead of the Ampère Glacier further supports glacier retreat at this stage (Frenot et al., 1997). The timing of the 1150 CE retreat at the Chamonix Glacier, combined with a previous, albeit shorter, retreat phase around 800 CE (Figures 10f, g), fits well with reconstructed patterns of Indo-Pacific sea-level pressure which suggests a two-stage Medieval Climate Anomaly in the SH, between 800–1100 CE and 1100–1300 CE, respectively (Goodwin et al., 2014). A change in the grain-size and carbon-to-nitrogen (C/N) ratio of central Chilean lake sediments (38°S) between 950 and 1250 CE (Kaiser et al., 2026), as well as elevated dust particle deposition in West Antarctica (79°S) around 1050–1400 CE are both suggested to represent a southward migration of the SHW (Koffman et al., 2014). Peaks in wind intensity in southern Patagonia (52°S) around 900 CE based on magnetite fluxes recorded in lake sediments support a poleward SHW displacement at this stage (Lisé-Pronovost et al., 2015). However, the calibrated age uncertainty at this interval in our records warrants caution when drawing parallels with other records (>260 years at 2σ). After 1200 CE, the Chamonix Glacier re-advanced, maintaining substantial ice extent until the late 16th century CE. At the same time, the Ampère Glacier advanced before following a similar pattern of variability with gradual retreat until present-day, although with a much higher amplitude of variability (Chassiot et al., 2024). A SHW reconstruction using hydrogen isotope ratios of lipids in lake sediments from South Georgia (54°S) is suggested to reflect windier conditions and negative SAM-like conditions the last millennium (van der Bilt et al., 2022). At the same time, another record of SHW changes, using hyperspectral ratios captured in lake sediments from Macquarie Island (54°S), another sub-Antarctic landmass located at the same latitude (Figure 1), interpret windier late Holocene conditions as indicative of a poleward shift of the SHW (Saunders et al., 2018), characteristic of positive SAM-like conditions. This is also recorded in conductivity and diatom data from another Macquarie Island lake, that show stronger winds between 1150 and 1850 CE (Deng et al., 2026).

5.5 Gradual decline around 1600 CE

From the late 16th century and onward, general retreat at the Chamonix Glacier parallels a shift towards more positive SAM-like phases, as consistently reconstructed across multiple SAM indices (Villalba et al., 2012; Abram et al., 2014; Dätwyler et al., 2018; King et al., 2023), supporting the leading view that more consistently positive SAM-like phases preceded industrial-era warming (Hessl et al., 2017). The same pattern is observed in an ice core from Law Dome in eastern Antarctica where, after 1600 CE, there is a shift towards lower interdecadal variability of the mean sea-level pressure, which is directly correlated to the SAM phasing (Goodwin et al., 2004). Furthermore, changes in the wind belt around this time interval is also documented in the SHW reconstruction from South Georgia (van der Bilt et al., 2022) as well as in a lake sediment conductivity-based SHW record from Marion Island (46°S) (Figure 1f) (Perren et al., 2020). Notably, greater glacier retreat is inferred after 1800 CE (Figure 10f), aligning with rising hemispheric temperatures (Neukom et al., 2014), a pattern also documented locally (Bellet et al., 2026). Nevertheless, the Chamonix Glacier had already experienced significant and long-term retreat before the advent of observational records in the 1960s CE (Favier et al., 2016). A smaller glacier re-advance can be deduced in the HEP-19–02 composite record that may be linked to the last period of glacial sediment input in Lake Aphrodite (Figure 10e), between 1790 and 1860 CE (Sleire et al., 2025). This coincides with the late Holocene MAAT minimum at the Kerguelen Islands, which occurred around 1840 CE (Bellet et al., 2026), yet the increase in clastic input is seemingly minimal compared to previous advances. It is possible that the decrease in MAAT was not sufficiently long to facilitate extensive glacier expansion. However, only one of our records cover this timeframe (Figures 10f, g) and, furthermore, we cannot exclude that chronological uncertainties explain the apparent offset.

6 Conclusion

Collectively, the glacier reconstructions of the Agassiz and Chamonix glaciers presented here alongside regional palaeoclimatic evidence, suggests that the Common Era climate around the Kerguelen Islands was characterised by extended periods of generally cold and wet (negative SAM-like) conditions that favoured glacier advances. However, short-term retreat phases as well as a longer, more gradual decline in glacier extent from around 1600 CE to the present is superimposed on this trend and may be linked to positive SAM-like anomalies that warmed and dried the atmosphere. The glacier sediment signal during the 600–750 CE event indicate extensive retreat of the Chamonix Glacier. Concurrently, large-scale remobilisation of the Chamonix Glacier coincided with negative SAM-like conditions. During these periods, the glacier extent was comparable to the largest extents documented in the available moraine records. The gradual decline in glacier activity around 1600 CE and until present-day supports the consensus that more positive SAM-like anomalies preceded industrial-era warming. If the prevailing positive SAM phasing persists under ongoing anthropogenic forcing (Deng et al., 2022) we, in line with previous studies, e.g., Verfaillie et al. (2021), anticipate continued, and potentially accelerated glacier disappearance on the Kerguelen Islands as well as other glaciated regions of the SH mid-latitudes.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. Generative AI was used solely to assist with the visual rendering of an author-created map. The underlying DEM, geographical information, and all scientific features were created and verified by the authors. No scientific data or geographical features were generated or modified using AI.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2026.1877993/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Joseph Licciardi, University of New Hampshire, United States

Reviewed by: Yafei Zou, Zhengzhou University, China

Federica Donda, National Institute of Oceanography and Applied Geophysics, Italy