Sample collection and laboratory analyses

Sediment cores extracted from Ardley Lake (ARD) and Yanou Lake (YAN) were analysed for the concentration of elements associated with changing inputs of penguin guano and volcanic ash deposits using complementary multi-proxy biogeochemical techniques described in this section and in Supplementary Methods.

Chronologies for the lake sediment cores were established using 18 (ARD) and 15 (YAN) AMS radiocarbon (14C) ages from, in order of preference: (1) moss macrofossil layers (consisting of hand-picked fine strands of the aquatic moss Drepanocladus longifolius (Mitt.) Paris, but also occasional layers of Campylium polygamum (Schimp.) Lange & C.E.O. Jensen, and some unidentifiable/mixed species moss fragments—considered more likely to have been reworked; (2) terrestrial and/or lacustrine algae; (3) other intact macrofossils and sub-fossils, including bones (bone-collagen, where extractable); (4) other (macro)fossil fragments; (5) organic-rich bulk sediments and, near the base of each core and as a last resort, (6) bulk glaciolacustrine or glaciomarine sediments (Supplementary Table 2). Bulk sediments were only dated where macrofossils were absent, while paired macrofossil or bone-collagen and bulk sediment samples were measured in the surface sediment and wherever present to check for any systematic offsets between ages obtained from different carbon sources (see Supplementary Methods, Supplementary Notes 5 and 6 for more details).

Measured radiocarbon ages from samples shown in Supplementary Table 2 were calibrated using SH13 and MARINE13 calibration curves and age-depth models were generated using Bayesian age-depth modelling techniques (Fig. 2, Supplementary Figs 3,9, Supplementary Methods and Supplementary Notes 5–7). All the ‘as measured’ (uncalibrated) radiocarbon age data shown in Supplementary Table 2 were used as input data for final age-depth model runs (ARD-M5, YAN-M4, ANVERS-M1, where model run number is indicated by the -M suffix). The weighted mean basal age of the ARD core was 8,750 cal a BP [8,410–9,230 min –max. 95% confidence age range]. Error analysis shows that whole record mean 95% confidence age-depth model uncertainties, rounded to the nearest 10 years with 95% confidence minimum to maximum age ranges (in years at depths in cm) shown in square brackets are: ARD-M5: 410 years [0.6 years at 0 cm–1,430 years at 335 cm], YAN-M4: 600 years [160 years at 7.7 cm–1,990 years at 340 cm], and ANVERS-M1: 1,420 years [590 years at 0 cm–1,810 years at 375 cm]. Equivalent values for the late Holocene in each record are: ARD-M5: 640 years [0.6 years at 0 cm–750 years at 20 cm], YAN-M4: 680 years [160 years at 7.7 cm–760 years at 2.7 cm] and ANVERS-M1: 1,370 years [590 years at 0 cm–1,600 years at 155 cm] (Supplementary Note 7).

Sub-samples for carbon, nitrogen, XRF and ICP-MS (quantitative, dry mass) bio-element analyses were taken at 1 cm intervals from ARD and YAN lake sediment cores, lyophilized and ground with an agate ball mill and manually with an agate pestle and mortar. ARD samples were analysed for total carbon (TC), total sulphur (TS) and total nitrogen (TN) using a CNS analyser (vario EL Cube, Elementar, Germany) equipped with a solid-state infrared and a heat conductivity detector. TC and TS measurements were conducted on YAN sediments by means of an ELTRA CS analyser. For data replication and comparison purposes, the total inorganic carbon (TIC) of 69 (out of 385) samples from both cores was determined coulometrically by a CM 5012 CO 2 coulometer coupled to a CM 5130 acidification module (UIC, USA) while total organic carbon was then calculated as the difference between TC and TIC (TOC%=TC%−TIC%). Owing to the high correlation between TC and TOC (R2>0.9997) and the negligible TIC concentrations (av. 0.10 and <0.01 mass%) compared to TC values (av. 6.2 and 2.5 mass%) in both cores, TC is considered to reflect the amount of organic carbon in both lake sediments. Quantitative XRF analysis for major and trace elements (SiO 2 , Al 2 O 3 , CaO, K 2 O, Na 2 O, P 2 O 5 , As, Ba, Cu, Co, Ni, Sr, Y, Zn, Zr) was carried out with a conventional wavelength dispersive X-ray fluorescence (WD-XRF) spectrometer (Philips PW 2400). For WD-XRF, glass beads were prepared following standard procedures74 and measurements undertaken in random order to avoid artificial trends. Trace element analysis (Cd, REE) of selected samples was performed by Inductively Coupled Plasma Mass Spectrometry (ICP-MS, Element 2 mass spectrometer, Thermo Scientific, Germany) at 2,500-fold dilution. For additional details concerning ICP-MS measuring conditions, see Schnetger75. Selenium was determined on acid digestions by graphite atomic absorption spectrometry (G-AAS) using a Unicam 939 QZ AA spectrometer and a Zeeman-effect background correction. A Milestone DMA-80 Direct Mercury Analyser was used for the measurement of mercury via cold vapour atomic absorption spectroscopy (CV-AAS). ICP-MS and G-AAS standard acid digestion procedures are described in Supplementary Methods.

Following established procedures, we used visual descriptions, smear slides and contiguous micro-XRF (μ-XRF) scanning, at 200 μm and 2 mm intervals, to determine the precise position of tephra deposits (see Supplementary Methods, Supplementary Note 9, Supplementary Figs 3,7,10,12 for details). Electron Probe Micro-Analysis (EPMA) of 165 glass-shards from the most prominent T4–T7 tephra layers in the YAN, ARD and Beak 1 Lake records (Beak Island; Layers T a–e (ref. 63)) was used to determine the eruption characteristics and source of mid-late Holocene eruptions that likely had the biggest environmental impacts. Results were compared with glass-shard analyses from age-equivalent tephra layers in marine cores PC460/461 from the Scotia Sea (new data; this study) and similarly analysed data in Moreton and Smellie23 from the Scotia Sea, northern Weddell Sea and Boyd Strait, and in Fretzdorff and Smellie22 from the Bransfield Basin (see Supplementary Methods, Figs 20,21, Supplementary Table 10).

For the YAN-GDGT reconstructed temperature record, temperature-sensitive glycerol dialkyl glycerol tetraethers (GDGTs) biomarkers were extracted from 41 samples using a microwave-assisted solvent extraction76,77. Contiguous 1 cm subsamples were taken in the top 20 cm of the YAN record, and at 2 cm intervals between 20 and 38 cm and 188 and 210 cm core depth (dated to c. 6.1 cal ka BP). Freeze dried and homogenized sediment samples weighing 0.2–4.3 g were microwave-extracted in DCM:Methanol (3:1, v/v). The total extracts were saponified and GDGTs isolated following ref. 76. Prior to analysis the GDGT extracts were filtered through a 0.2 μm Whatman PTFE filter. GDGT analysis was undertaken using an Acquity Xevo TQ-S (triple quadrupole with step wave; Waters Ltd.) LC-MS set up with an atmospheric pressure chemical ionization source (Ion saber II) operated in positive ion mode. Analytical separation was achieved using a Grace Prevail Cyano HPLC column (3 μm, 150 × 2.1 mm i.d.) fitted with an in-line filter (Waters Acquity UPLC in-line filter, 0.2 μm) at 40 °C using a binary solvent gradient where eluent A was hexane and eluent B was propanol. The flow rate of the mobile phase was 0.2 ml per minute with a gradient profile of 99% A 1% B (0–50 min); 98.2% A 1.8% B (50–55 min); 90% A 10% B (55–65 min) and finally 99% A 1% B (66–80 min). The LC-MS settings were: source offset 50 V, capillary 1.5 kV, desolvation temperature 200 °C, cone voltage 30 V, desolvation gas (N 2 ). Detection was achieved using selected ion monitoring of targeted [M+H]+ ions (dwell time 50 ms). The target ions were m/z 1,302, 1,300, 1,298, 1,296 and 1,292 for the isoprenoid GDGTs (isoGDGTs) compounds and 1,050, 1,048, 1,046, 1,036, 1,034, 1,032, 1,022, 1,020 and 1,018 for the branched GDGT (brGDGTs) compounds. GDGTs were identified and integrated using MassLynx software (v.4.1) and GDGT-derived temperatures were calculated using the Antarctic and sub-Antarctic GDGT-temperature calibration77 (see Data Analysis section). The Antarctic and sub-Antarctic GDGT-MSAT surface calibration data set (comprising 32 sites in total) includes surface sediments from Yanou Lake, Ardley Lake and three other lakes from Fildes Peninsula, along with two further lakes from Potter Peninsula and four lakes from the Trinity Peninsula, north-eastern AP (including Beak Island).

Diatom analysis of the ARD and YAN records is as described in ref. 45 and summarized in Supplementary Fig. 14, while Anvers Shelf marine core (ANVERS) diatom analysis and chronology are described in Supplementary Methods and summarized in Supplementary Fig. 15.

Data analysis

Bio-element assemblages are immobile in (Antarctic) lake sediments43 and guano input to sediments leads to the formation and preservation of stable phosphates, such as struvite (Mg(NH 4 )PO 4 x 6 H 2 O), leukoposphite and, in particular, hydroxylapatite (Ca 5 (PO 4 ) 3 (OH)), which is one of the dominant compounds found in ornithogenic soils/sediments on King George Island. During the precipitation of apatites, an exchange between Ca2+, PO 4 3−, F− and OH− and elements such as Ag, Br, Ba, Cd, Cu, Cr, I, Na, Mg, Mo, Pb, S, Se, Sr, U, V, Y and Zn is possible, coupled to microbial-mediated degradation of solid phases78 (see Supplementary Note 1 for more details). These trace elements, which have naturally higher concentrations in penguin guano, are enriched in these phosphate phases and then immobilized in soils and sediment profiles during this substitution process. While the high correlation between Mg and P in ornithogenic soils from Vestfold Hills has been used to indicate the presence of struvite79, their lack of correlation in the Ardley Lake sediments means that Mg is likely derived from the bedrock lithology (Supplementary Fig. 3, Supplementary Table 4), and hydroxylapatite (Ca 5 (PO 4 ) 3 (OH)) formation reflects guano input (Supplementary Fig. 6).

The average relative proportion of guano or ornithogenically derived sediment in the Ardley Lake sediment matrix (F o.sed. ) was estimated by using the mixing equation (equation (1)), modified after Shultz and Calder80:

(El/Al) smp is the element/aluminium ratio of selected bio-elements (Cu, P, Sr, Zn) in Ardley Lake sediments and (El/Al) bgd (mean Ardley Island bedrock51) and (El/Al) o.sed. (guano-bearing ornithogenic sediment and soils49) represent the respective ratios in both end members. In order to minimize misinterpretation of this proxy, we calculated the average fraction of ornithogenic sediment and soils in the lake sediments using a combination of four chemically unrelated El/Al ratios. This provided a buffering effect on possible variations in ornithogenic soil and bedrock composition over time. In the resultant F o.sed. percentage plotted in Fig. 4b, for example, F o.sed. =20%, means that 20% of the sediment sample is eroded ornithogenic sediment-soil and 80% eroded bedrock. Theoretically, all F o.sed. values >0% provide evidence for guano input above the local bedrock-terrigenous background level, but we used a weighted mean F o.sed. value of >10% at the cut-off for CONISS-defined geochemical zone to indicate the sustained presence of penguins around Ardley Lake (Table 1, Supplementary Table 7). We calculated post-eruption colony recovery intervals as the time taken until the percentage F o.sed. value first returned to >10% within the subsequent guano phase to represent the return of a sustained penguin presence around Ardley Lake. Measurement errors of <1% mean our definitions provide very conservative buffers of penguin presence during guano phases 1–5 (GP-1 to GP-5) (Table 2).

Penguin population modelling results shown in Fig. 5 utilized dry mass F o.sed. . accumulation rates (F o.sed. DMAR in Fig. 4b), which were calculated following standard procedures81 (that is, multiplying F o.sed. percentages by the dry mass accumulation rate, calculated from sedimentation rate and dry mass density data). We also used published penguin population parameters40,41 and the ‘most-realistic’ sediment focussed ellipsoid basin accumulation models81, assuming a <30 ° slope-angle accumulation area within Ardley Lake of 5,682 m2 (out of a maximum 7,274 m2 lake basin area) (see Supplementary Fig. 2). Further accumulation area scenarios, equations and references can be found in Supplementary Note 8. We considered the shape of Ardley Lake basin to approximate between a steep-sided hyperboloid basin and an ellipsoid basin (as below).

To summarize, the basin-wide guano influx rate, S G , at time t, and total Ardley basin maximum core depths of z (water depth + sediment depth, in cm) at time t, was calculated as follows:

For an ellipsoid-basin:

and

For a hyperboloid basin shape:

where z=core depth at time t; z m =maximum depth of the original basin (that is, current lake depth + sediment record depth); =ratio of the mean depth to the maximum depth of the original basin (z m ); =guano dry mass accumulation rate=F o.sed DMAR in g cm−2 a−1.

Assuming Ardley Lake best approximates to an ellipsoid basin shape (that is, the ‘more likely’ basin shape scenario), the total basin-wide guano accumulated (I G ) over the ‘more likely’ <30 degree slope accumulation area (A <30 ) within Ardley Lake (Supplementary Fig. 2), corrected for the effects of sediment focussing, is then:

where ; = the value for each 100-year interval obtained from equation (2).

Using guano production rates of 84.5±21.1 g (25% error estimate applied) per penguin per day on Ardley Island, and a mean density of 0.31±0.19 gentoo penguins per m2 (for Signy Island)40,41, we estimated the total amount of guano delivered by erosion from the catchment into Ardley Lake as the guano yield (G y ) and the area required for the reconstructed population. As the main breeding season and ice/snow-free periods, when active erosion from catchments to lakes occurs on Ardley Island, last for approximately 3–4 months each year (91.3±22.8 days; 25% error estimate; equation (6) and Supplementary Table 9), we calculated the maximum amount of guano deposited per penguin (p−1), per year (a−1) inside the catchment and delivered to the lake as the maximum possible guano yield (G y-100% ) as follows:

where =mean occupation time per year±estimated error=365.25/4 or 91.3±22.8 days (25% error applied); =amount, in grams (g) of guano deposited per year (a−1) per penguin (p−1).

Several studies have shown that catchment erosion and deposition in small closed-basins with a single catchment and non-complex inflow characteristics (for example, Ardley Lake) can be approximated by a linear relationship (see Supplementary Note 8 for references); thus, a reasonable approximation of the total amount of sediment or, in this case, guano delivered (G D ), by erosion (G e ), into the Ardley Lake basin, can be obtained, simply from ratio of lake-area to catchment-area, expressed as an estimated percentage:

where G D is the ratio, expressed as a percentage, of the total guano produced in the catchment (G y-100% ) and the total guano-eroded from the catchment (G e ) into the lake, which, here we approximate to the lake area: lake catchment ratio, .

For Ardley Lake, the amount of guano delivered in grams (g) per year (a−1) per penguin (p−1) in the 11% guano (G y–11% ) yield scenario is given by:

However, small lakes with restricted inflow and/or low rates of catchment erosion or slow sedimentation rates (that is, most typically high latitude, frozen lakes, for example, Yanou Lake during the late Holocene) are better represented by an exponential sediment delivery model, (ref. 81), which, for Ardley Lake, produces as follows:

This allowed us to estimate the number of penguins (P) present in the catchment at time (t) for three guano-yield scenarios, :

as follows

Since the catchment colony would have produced and/or deposited guano outside the catchment area (for example, on foraging trips), and not all guano deposited in catchment areas is currently eroded into lake due to various factors, for example, utilization by vegetation, we consider the G y–100% to be very unlikely. Therefore, we consider the G y–11% and to best represent a possible ‘most likely’ minimum population scenario as follows:

As a reality check, the G y-11% and scenario recreated the near-absence of a penguin colony (<10 penguins) in the modern day catchment. Conversely, we consider that the G y-0.96% and scenario represents the ‘most likely’ maximum population scenario:

These two end-member scenarios are shown in Fig. 5, but since sedimentation rates in Ardley Lake have varied significantly through time, it is possible that, when erosion rates are high, for example, between c. 5.5 and 4.5 cal ka BP, the Ardley Island colony population would be better estimated by a guano-yield scenario somewhere between 11 and 100%. To further validate our findings, we also calculated the catchment area required by the reconstructed colony population using published gentoo penguin density values40 and calculated mixed-species density values, partitioned according to published modern-day Ardley Island gentoo-Adélie-Chinstrap species compositions41 and the corresponding species density values40 (Supplementary Fig. 19b, Supplementary Table 9).

The ‘most-likely’ minimum guano-delivery (yield) scenario assumes that 11% of the guano produced by penguins in the Ardley Lake catchment was eroded into Ardley Lake. This scenario produced a peak guano phase GP-4 colony population of 1,024±574 penguins, which occupied an area of 1,679–3,303 m2 out of the 66,249 m2 Ardley Lake catchment area (based on partitioned mixed species minimum density–gentoo penguin maximum density values40). Meanwhile, the 0.96% guano-yield scenario produced a maximum peak GP-4 colony population of 11,723±5,104, which corresponds to an occupied area of 20,686–37,815 m2 (Supplementary Fig. 19). Therefore, the c. 66,250 m2 Ardley Lake catchment area is easily large enough to accommodate all modelled population scenarios shown in Fig. 5. Assuming a hyperboloid basin shape for Ardley Lake would result in an average 33±25% reduction in the reconstructed population across all guano delivery scenarios.

It is important to note that the apparent increase in the amount of guano deposited into Ardley Lake over the last 500 years is less than the significant 10% F o.sed. level we used to determine colony presence elsewhere in its Holocene record, and only equates to a total of 65±48 penguins (using the 11% guano-delivery scenario). This also illustrates why the Ardley Lake guano phase colony was too small after c. 2,500 years for us to determine whether volcanic activity had a significant impact.

Correlation analysis of geochemical data was undertaken using R 2.15.2 (R Foundation for Statistical Computing). Hierarchical R-mode cluster analyses were conducted with the R package ‘Pvclust’ (version 1.2-2) using average linkage and a correlation-based dissimilarity matrix (see Supplementary Methods for software references). Based on multi-scale bootstrap resampling (number of bootstraps: 10,000), approximately unbiased P values were further calculated to assess the uncertainties in cluster analysis. We used constrained incremental sum of squares (CONISS) stratigraphically constrained cluster analysis with broken stick analysis to define Diatom Zones (D) in the ARD, YAN and ANVERS records, Geochemical Zones (GZ) in the ARD and YAN records, and Guano phase GP-1 to GP-5 boundaries in the ARD record. Downcore cluster-zonation was undertaken on ecologically grouped diatom percentage data and square-root transformed geochemical data using the R packages vegan and rioja (see Supplementary Methods for software references). Downcore trends in diatom and geochemical data sets using Principal Components Analysis (PCA) were undertaken using square-root transformed percentage XRF data and μ-XRF (Total Scatter Normalised (TSN)) percentage geochemical data and percentage diatom count data.

For climate and sea-ice comparisons, the Yanou Lake reconstructed GDGT temperature (YAN-GDGT) data (grey plot, with RMSE errors shaded grey in Supplementary Fig. 16a) and the Maxwell Bay measured TOC data (sea-ice proxy record)52 (grey plot in Fig. 4c) were smoothed using polynomial negative exponential regression analysis. Results were similar to smoothing by LOESS first-order polynomial regression with tri-cube weighting (not shown). Reconstructed GDGT temperatures reflect MSAT (December, January, February or DJF) in the shallow and well-mixed open lake water environment of Yanou Lake76,77. For the YAN record, whole YAN-GDGT data set 5th and 95th percentile outliers (<−0.33 °C and >8.61 °C) 11.34 °C and 13.84 °C at 1,164 cal a BP and 3,135 cal a BP, respectively, were excluded by smoothing analysis (Fig. 4a,c, Supplementary Fig. 17a). These values were also greater than the 10.3 °C maximum limit of the surface sediment-MSAT Antarctic and sub-Antarctic GDGT-MSAT calibration (ANT-GDGT) data set77. All reconstructed GDGT temperatures were greater than the −2.2 °C minimum limit of the ANT-GDGT calibration data set. The YAN-GDGT <10 °C data set was used in 6–0 ka weighted mean temperature anomaly calculations, statistical analyses and hypothesis testing. The mean±1σ 6–0 ka reconstructed YAN-GDGT <10 °C data set weighted mean MSAT value of +2.39±2.67 °C [95% confidence interval: 1.52–3.25 °C] (Bootstrapping, n=10,000, +2.38±2.61 °C [1.52–3.25 °C]; 5th/95th percentile=−0.51/+7.96 °C; weighted by Antarctic GDGT calibration data set RMSE value of 1.45 °C; Supplementary Table 7) encompasses the modern-day (1968–2015 CE) observational mean summer (DJF) temperature from Bellingshausen Research Station on Fildes Peninsula of +1.09±0.56 °C (2σ 95% measured range=−0.03 to 2.21 °C; n=48, 1968–2015 CE, where year refers to December) (Marshall, pers. comm.; SCAR-READER database: http://www.antarctica.ac.uk/met/READER/). On a regional to global-scale modern-day reconstructed lake sediment GDGT-temperatures broadly reflect changes in air temperature77. At the local, small lake scale, water temperature might decouple from observed MSAT due to factors such as the thickness and length of seasonal lake ice-cover and the volume of meltwater input, both of which suppress lake water temperatures. Conversely, since water bodies store excess heat generated during extended ice-free seasons, small and/or shallow lakes (<10 m deep) reach temperatures of 6 °C or more (for example, Sombre Lake, Signy Island)82. As the focus of this paper is penguin colony reconstruction, we investigate these aspects further in a forthcoming paper.

After smoothing, all climate (YAN-GDGT, JRI29 and PD48); and sea-ice (Anvers Shelf and Maxwell Bay52) data sets were resampled at 100-year intervals (see Supplementary Methods). We then used >6–0 ka upper bound on the weighted mean (95% confidence level) temperature, temperature anomaly, diatom open water ratio and smoothed TOC- and diatom-based sea-ice proxy data (Fig. 4) to define ‘warm’ and ‘sea-ice free’ (open-water) periods. We chose the 6–0 ka period and values greater than the upper bound on the weighted mean at 95% confidence level of the weighted mean as this provides the best longer-term assessment of the mean temperature and variability that penguin colonies on Ardley Island would have encountered. It also corresponds to the time period covered by the terrestrial Yanou Lake GDGT-MSAT record. Using the Pre-industrial era (1,000–250 cal a BP) and its upper bound on the weighted mean value (95% confidence) produced similar timings for ‘warm’ and ‘open-water’ intervals.

To compare variations in cross-Peninsula climate during guano and non-guano phases, we undertook normality tests and produced a set of descriptive statistics for the whole Ardley F o.sed. data set (weighted by 1σ error), the YAN-GDGT <10 °C temperature data set, the JRI temperature anomaly data set29, and the PD 0–200 m sea-surface temperature (PD-SST) record48 (weighted by published errors, where available). Using α=5% and 10% significance levels for three hypothesis test scenarios A–C shown in Supplementary Table 8, we first performed Fisher’s F-test Variance Ratio analysis. Since some data sets were not normally distributed, we then used the Mann–Whitney U-statistic to test the hypothesis that temperature and temperature anomaly distributions in the YAN-GDGT <10 °C, JRI, PD-SST records were statistically different during the five guano-influenced phases (Scenario 1A in Supplementary Table 8) compared to the eight non-guano phases (Scenario 2A). We then repeated this process comparing guano phases of the mid-Holocene (8.2–4.2 ka; 1B) to those of the late Holocene (4.2–0 ka; 2B), and for late Holocene guano phases (1C) versus non-guano phases (2C). We also undertook similar hypothesis testing using the Maxwell Bay sea-ice and Anvers Shelf open-water data sets to determine whether sea-ice conditions within guano and non-guano were statistically different in the three scenarios (A–C) described above. Results are summarized in Supplementary Table 8.

Apart from one recent minor eruption, all visible (that is, >2 mm thick) Holocene-age tephra deposits in lake, marine and ice core records of the northern Antarctic Peninsula (N-AP), SSI and Scotia Sea have been linked to VEI=3 or VEI>3 eruptions from Deception Island20,21,22,23,27,28,29,58,63,83,84 (Supplementary Note 9, Supplementary Figs 20, 21, Supplementary Tables 3, 10). Tephra from even relatively small 1967 and 1970 CE Deception Island eruptions (VEI=3) reached more than 150 km (ref. 22) and is present in the JRI ice core29. Neither the 1976 or 1970 eruptions formed visible ash layers near the top of the ARD core, probably because only up to 1 mm of ash was deposited on Fildes Peninsula21. Therefore, we broadly equated major VEI=3 eruptions from Deception to visible tephra deposits <1 cm thick in the YAN and ARD records, and airfall ash deposits >1 cm thick to VEI>3 eruptions. Since the post-deglaciation (after c. 8–7.5 cal ka BP) lithogenic deposition style in YAN and ARD is predominantly fine-grained silt-clay, the sand dry mass accumulation rate (DMAR) provided a useful proxy for eruption size and the volume of tephra deposited onto Ardley Island and King George Island during VEI>3 eruptions (Fig. 4a,b, Supplementary Figs 12,20,21), while elevated Ca/Ti ratios and PCA2 values from combined 2 mm and 200 μm μ-XRF scanning data and the reliable single-species aquatic moss age-depth model provided precise stratigraphic constraints on the positions (and ages) of tephra deposits in the YAN record (Supplementary Note 9, Supplementary Figs 9–11,16, Supplementary Table 3).

Data availability

Key data sets for this study are in Supplementary Data files 1–7. All data produced by this study and associated IMCOAST/IMCONET projects are available from the corresponding author (sjro@bas.ac.uk) and the NERC Polar Data Centre (BAS) (https://www.bas.ac.uk/data/uk-pdc/) and Pangaea Data Publisher website (AWI) (https://www.pangaea.de/). Antarctic location maps and spatial data in Fig. 1 are the LIMA (Landsat Image Mosaic of Antarctica) downloaded from https://lima.usgs.gov/index.php and the SCAR Antarctic Digital Database http://www.add.scar.org/, used under a Creative Commons license. The basemap in Fig. 1a and Supplementary Figs 20,23 is the ESRI Ocean Basemap (http://goto.arcgisonline.com/maps/Ocean_Basemap) used under a BAS-NERC software license. Satellite images and image derivatives in Figs 1c,d and 6 and Supplementary Figs 2,22,23 are copyright of DigitalGlobe, Inc. All Rights Reserved (Catalogue numbers: 1030010020C0C900 and 11MAR04115702), reproduced at low resolution in a static format under a BAS-NERC educational license, and are not included in the Creative Commons license for the article.