Antarctic ice sheet discharge driven by atmosphere-ocean feedbacks at the Last Glacial Termination

Reconstructing the dynamic response of the Antarctic ice sheets to warming during the Last Glacial Termination (LGT; 18,000–11,650 yrs ago) allows us to disentangle ice-climate feedbacks that are key to improving future projections. Whilst the sequence of events during this period is reasonably well-known, relatively poor chronological control has precluded precise alignment of ice, atmospheric and marine records, making it difficult to assess relationships between Antarctic ice-sheet (AIS) dynamics, climate change and sea level. Here we present results from a highly-resolved ‘horizontal ice core’ from the Weddell Sea Embayment, which records millennial-scale AIS dynamics across this extensive region. Counterintuitively, we find AIS mass-loss across the full duration of the Antarctic Cold Reversal (ACR; 14,600–12,700 yrs ago), with stabilisation during the subsequent millennia of atmospheric warming. Earth-system and ice-sheet modelling suggests these contrasting trends were likely Antarctic-wide, sustained by feedbacks amplified by the delivery of Circumpolar Deep Water onto the continental shelf. Given the anti-phase relationship between inter-hemispheric climate trends across the LGT our findings demonstrate that Southern Ocean-AIS feedbacks were controlled by global atmospheric teleconnections. With increasing stratification of the Southern Ocean and intensification of mid-latitude westerly winds today, such teleconnections could amplify AIS mass loss and accelerate global sea-level rise.

Melted ice samples retrieved from below the surface with a Kovacs 9cm diameter ice corer were 79 centrifuged at 2500 rpm, and the remaining particulate material was mounted on glass slides for 80 electron microprobe analysis. The slides were ground using silica carbide paper and polished using 81 9, 6 and 1 µm diamond suspensions to produce fresh sections of the glass shards. Electron 82 microprobe analysis was undertaken at the Tephrochronology Analytical Unit at the University of 83 Edinburgh. Single-grain analyses of ten oxides were performed on a Cameca SX-100 instrument 84 with 5 wavelength dispersive spectrometers. A 5 µm beam diameter was used and operating 85 conditions followed those in 5 . Lipari Obsidian and BCR2G were analysed as secondary standards. 86 A summary of the geochemical results is provided in Table S1. 87 Geochemical results reveal a bimodal composition for 282 m (Table S1) and single tightly-89 clustered trachytic populations for PH279 m and PH190 m (trachyte: PH282-a and rhyolite PH282-90 b) (Table S1). PH282-b, PH279 and PH190 show geochemical affinity to Mt Berlin products and 91 to tephras identified in the Siple Dome core 6 , Mt Moulton 7 and EPICA 8 . Similarity coefficient 92 (SC) analysis 9 for PH190, highlight strong matches to a number of Siple Dome and Mt Moulton 93 tephras ranging in age between 37 and 18 ka (Table S2). Several potential matches are revealed 94 by SC values ≥0.95, with one sample, SDMA-5951c set 1 exhibiting a SC of 0.99 (Table S2A). 95 Bivariate plots of samples with SC ≥ 0.97 for PH190 are shown in Figure S1. Although the data-96 sets overlap on several major elements we rule out SDMA-5635c on the basis of higher FeO tot 97 values, SDMA-9063b on higher SiO 2 values and slightly lower K 2 O values and SDMA-5785c on 98 lower SiO 2 values and slightly higher FeO tot values ( Figure S1). The data-set for SDMA-5951c set 99 1 exhibits a slightly wider range than PH190 but we believe that this is the most likely match. 100 SDMA-5951c set 1 is dated to 36.4 ka yrs in the Siple Dome ice core 6 . 101

102
The trachytic tephra PH279 most probably relates to volcanism in Marie Byrd Land 8 . Five tephras 103 have SC values ≥ 0.95 with PH279 (Table S2B) and the compositional data for these are plotted 104 on Figure S2. SDMA-9063b can be ruled out on the basis of higher SiO 2 values and EDC434.7 105 exhibits higher Al 2 O 3 compositions. Although the remaining tephras overlap on many of the major 106 elements, shards from PH279 are distinctly higher in Na 2 O than SDMA-5554c and SDMA-5694. 107 Although only summary statistics are available for WCM93-25, this is the only deposit to show 108 partial overlap with PH279 on Al 2 O 3 vs Na 2 O and SiO 2 vs Na 2 O plots ( Figure S2). On this basis, 109 we correlate PH279 with WCM-93-25 dated to 18.2 ka in the Mt Moulton record 7 . 110 The rhyolitic component of the tephra at 282 m (PH282-a) is tightly clustered with average SiO 2 112 values of 74.6%, FeO tot values of 2.1% and total alkali values of 9.7% (Table S1). We are unable 113 to suggest a potential correlation for this tephra. In fact very few tephras of rhyolitic compositions 114 are reported from the Antarctic ice-core records 6,8,10 , and we suggest this tephra may originate 115 from a volcanic source outside Antarctica. The trachytic population of PH282-b again shows 116 considerable scatter and probably relates to volcanism in Marie Byrd Land 8 , but we are unable to 117 pinpoint a tephra of similar age and composition. As outlined we suggest that PH282-a and -b 118 may relate to two separate eruptions in different volcanic systems. We suggest that PH282a and b 119 most likely relates to widespread volcanic events dated to ~17.8 ka years and found in Byrd 11 and 120 the WAIS divide 12 cores, which to date has not been geochemically typed but used as an Antarctic Further chronological control across the profile is provided by a comprehensive suite of trace 125 gas samples -carbon dioxide (CO 2 ), methane (CH 4 ) and nitrous oxides (N 2 O) -taken from depth 126 (>3 m) along the BIA transect. The trace gases were extracted and measured at CSIRO's Ocean 127 and Atmosphere facility in Melbourne, and aligned to published values reported from EPICA 128 Dome C, providing a range of possible age solutions, that together with the absolute constraint 129 provided by the tephra horizons, allows the development of a robust chronological framework that 130 can be tied directly to the isotopic profile through high-resolution GPR survey 2,3 . The available 131 constraints indicate the complete 800 m long Patriot Hills BIA transect spans ~50 to ~2.3 ka (main 132 text Figure 2). Trace gases were extracted from the ice, was undertaken at CSIRO Ice Lab using a 133 dry extraction 'cheese grater' and cryogenic trapping technique 13 , with minor alterations 14 . The trapped air samples were analysed by gas chromatography (GC) and the trace gas concentrations 135 are reported on the calibration scales maintained by CSIRO GASLAB 15 . 136

137
The concept behind the dating routine is that CH 4 is the most reliable gas, as it is hard to be 138 produced in-situ and it is generally not altered by post coring melting as much as CO 2 and N 2 O. 139 Therefore, we use "CH 4 only" to derive a first estimate of the age and age distribution. Then we 140 use CO 2 and N 2 O to better constrain the age distribution, only if they agree with the result provided 141 by "CH 4 only". The criteria used to decide whether CO 2 and N 2 O agree with the result from "CH 4 142 only" is that the ages estimated with the different gases differ either by less than 1000 years, or by 143 less than 5000 years and the CO concentration is less than 120 ppb (high CO concentrations 144 suggest possible in-situ production due to reaction of organics with hydrogen peroxide (H 2 O 2 )). 145 Furthermore to assess possible modern atmospheric contamination we undertook analysis for the 146 contaminant sulfur hexafluoride (SF 6 ) on a sample sub-set. The average concentration of eight 147 samples analysed for SF6 was about 5% of modern day atmospheric concentrations, and less than 148 2% for the two samples selected to develop the chronology (with values ranging from 0.01-0.14 149 ppt). These thresholds have been chosen based on the minimum age difference that we aim to 150 resolve and on the expectation of CO atmospheric values derived from past records 14 . 151

152
The chronological framework for the Patriot Hills BIA profile was defined from the trace gases 153 concentrations following nine key steps: 154 Step 1. Assign a value and an uncertainty to each sample: This is a critical step because the 155 width of the distribution of gas concentration values determines the range of possible ages 156 attributed to each sample.
The measured values have been corrected for systematic effects associated with the extraction 158 procedure (named blank correction) according to the numbers in the headers of Table S3. Two 159 different blank corrections have been used based on the period when the samples were extracted 160 or analysed (2013 or 2015) due to changes in the extraction procedure. The value attributed to 161 samples analysed in replicates is the average of multiple measurements. The uncertainty attributed 162 to them is the difference between measured values. The uncertainty attributed to samples that were 163 not measured in replicates is the average of the uncertainties shown in Table S3 (CH 4 = 27 ppb, 164 CO 2 = 4 ppm; N 2 O=2 ppb, see Table 4). Depending on the relative position and the slope of the 165 spline fit to CH 4 measurements, the samples are assigned to a specific period (Glacial, Transition, 166 Early or Late Holocene). These periods are converted into a start and an end of the possible gas 167 age distribution based on the presence of tephra though the BIA profile, as shown in Table S4. 168 While it was not possible to measure all samples in replicates, the agreement among age 169 distributions suggested by multiple gases provide great confidence in our attributed ages 170 Step 2. Produce a Gaussian distribution for each gas species: With mean value equals the 171 measured value and 1 sigma equals the uncertainty attributed (as shown in Table S4). This step 172 allows us to attribute a probability to each possible range of values because the probability of 173 obtaining a value in a certain concentration range is proportional to the area under the Gaussian 174 distribution. The Gaussian distribution is divided in 100,000 intervals and covers a 3-sigma range 175 (99.73 % of all possible values around the measured value). 176 Step 3. Fit the CH 4 , CO 2 and N 2 O records covering the last 120,000 years with a cubic spline: 177 The gas records are a compilation of records from different ice cores. For CH 4 , the EDC record 178 has been taken. For CO 2 , a compilation of records (mostly based on the EDC record 16 ) has been taken. For N 2 O, the Talos Dome record has been taken 17 . The cubic spline has been produced by 180 using the Matlab spline tool and choosing a time step of 10 years. 181 Step 4. Find the match between the Gaussian distribution and the spline fit for each species: 182 For each gas age interval (x values on the spline fit plot), our dating routine checks whether there 183 is a value of the Gaussian distribution that falls into the range of values (y values on the spline fit 184 plot) defined by the start and end of the interval. 185 If a match is found, the probability for that value from the Gaussian distribution is attributed to the 186 corresponding time interval. When more than one match is found for a specific gas age interval, 187 the sum of the probabilities is taken. This procedure provides a way of projecting the probability 188 from the Gaussian distribution onto the time axis ( Figure S3). In other words we convert a 189 "probability vs gas concentration" plot into a "probability vs gas age" plot for each gas species 190 which quantifies the likelihood of a sample being x years old for each measured gas species. To 191 represent a likelihood, the "probability vs gas age" plot (derived from each gas species) is 192 normalised to a total area of 1. 193 Step 5: Find the intersection with the spline fit to the atmospheric records: Each measured 194 value is projected to find the intersection with the atmospheric CH 4 spline fit. The most likely age 195 is the intersecting value with the highest probability (mean value of the CH 4 age distribution), 196 within the limit of the period that the sample falls in (Glacial, Deglacial, Early or Late Holocene, 197 defined by the absolute chronological information define by available tephras). The same is 198 performed with CO 2 and N 2 O separately. If the most likely CO 2 age agrees with the most likely 199 CH 4 age (according to the criteria described above), then a combined CH 4 /CO 2 age is derived, 200 based on the combined CH 4 /CO 2 age distribution. If the most likely N 2 O age agrees with the 201 combined CH 4 /CO 2 age, a combined CH 4 /CO 2 /N 2 O age is derived Step 6. Force order in the age sequence: The sequence of samples is ordered chronologically 203 from the oldest (PH-10) to the youngest (PH-800), meaning that e.g.: PH-270 must be younger 204 than PH-200 and older than PH-290 and so on. This is used to constrain the most likely age 205 attributed to each sample following the logical order based upon high-resolution GPR analysis 3 206 and the ages of the available tephra in the BIA profile. 207 Step 7: calculate the 1-sigma age range: Gas age ranges (68 % corresponding to 1 sigma) are 208 calculated by integrating the area corresponding to 68 % of the total area of possible gas ages 209 around the most likely age. 210 Step 8. Combine the likelihood from all species in one: The three gas species are assumed to 211 provide independent evidence of likelihood of gas age and can be multiplied to get a combined 212 likelihood (as when calculating the joint probability of independent events). We end up with a 213 combined "probability vs gas age" plot. The area under the curve quantifies how likely it is that a 214 sample is between x and x+∆x years old. Finally, the dating routine is run for all samples 215 independently. This is summarised in Figure S4, the full sequence of samples, and can be seen 216 against the EDC record of trace gases in Figure S5 18 . 217 218

Constructing a geochronological framework across the profile 219
The final age model between the unconformities D1 and D2 can be seen in Table S5. The isotopic 220 profile was corrected using the modelled bubble close off or ∆age from the WDC core, which has 221 the most tightly constrained ∆age of any other regional ice core record 19 . The six well constrained 222 tie-points that fall between the unconformities D1 and D2, lie across the Last Glacial Termination, 223 and in conjunction with the age of well constrained volcanic horizons, recorded in other Antarctic 224 ice cores (see section 2A) allow us to construct a robust chronology for the Patriot Hills BIA, and particularly between unconformities D1 (247 m) and D2 (360m), which spans the period ~23 ka 226 to 11 ka years (main text Figure 2 and Figures S5 and S6). Ocean Water 2 (VSMOW2).

Earth system and ice sheet modelling 271
Transient simulations of the last deglaciation was performed with the Earth System model 272

Discussion of dynamic behaviour of WAIS across the Weddell Sea
The dynamic nature of the surface profile changes suggested by our data and predicted by ice-317 sheet model simulations reflects sensitivity of the AIS to ocean forcing (main text Figure 3). This 318 together with evolution of the basal substrate and hydrology, allows ice streams -the narrow 319 conduits of fast flow which control the mass balance of the Antarctic ice sheet -to switch on and 320 off. The ability to capture this evolving dynamic ice stream behaviour is critical, and enables 321 detailed reconstruction of periods of rapid drawdown and flow direction changes to be achieved 322 ( Figure S8 The question over ice advection from outside Horseshoe Valley is crucial, as it impacts the 336 interpretation of the elevation changes recorded in the ice exposed at the Patriot Hills BIA, and 337 therefore any WAIS ice dynamic changes predicted. Key to understanding this are insights into 338 regional ice sheet dynamics derived from the interpretation of airborne radar, specifically radio shown (see Table S2). shown (see Table 2B