Pliocene deglacial event timelines and the biogeochemical response offshore Wilkes Subglacial Basin, East Antarctica

a Dept. Earth Science and Engineering, Imperial College London, UK b Grantham Institute – Climate Change and the Environment, Imperial College London, UK c Antarctic Research Centre, Victoria University of Wellington, Wellington, New Zealand d Binghamton University, New York, USA e Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Japan f Instituto Andaluz de Ciencias de la Tierra, CSIC-Universidad de Granada, Spain g Dept. of Geology, University of Otago, Dunedin, New Zealand h Dept. of Marine Science, University of Otago, Dunedin, New Zealand


Introduction
Ice grounded below sea level in the East Antarctic Ice Sheet (EAIS) has a potential sea level equivalent of ∼19 m (Fretwell et al., 2013). Such ice is predominantly contained within the Aurora Subglacial Basin, the Recovery Basin, and the Wilkes Subglacial Basin (Fig. 1a). Modern observations in the vicinity of the Aurora Subglacial Basin document significant glacier retreat (Miles et al., 2016) and basal melting driven by ocean warming (Rintoul et al., 2016), while recent modeling work suggests that the Recov-  Cook et al., 2013). Topographic map of the Wilkes Subglacial Basin (from Bedmap2; Fretwell et al., 2013) with simplified lithologies and their isotopic characteristics. Areas of inferred Jurassic Ferrar Large Igneous Province (FLIP) from airborne geophysics are shown in hatched marking (Ferraccioli et al., 2009;Frederick et al., 2016;Studinger et al., 2004). CB denotes the Central Basin. (For interpretation of the colors in the figure, the reader is referred to the web version of this article.) Previous studies within the vicinity of the Wilkes Subglacial Basin have utilized detrital sediment provenance, grain size, icerafted debris accumulation, grounding line advances and retreats, and sea ice and temperature reconstructions to reveal a dynamic picture of the early to middle Pliocene and middle Miocene ice margin, with substantial ice retreat into the basin during warmer times, and ice advance onto the outer shelf at colder times (Armbrecht et al., 2018;Cook et al., 2013Cook et al., , 2017Orejola et al., 2014;Patterson et al., 2014;Pierce et al., 2017;Reinardy et al., 2015;Sangiorgi et al., 2018). Substantial ice retreat into the Wilkes Subglacial Basin during the Pliocene and Miocene is also indicated by models that incorporate complex flow regimes of ice streams, capture marine ice sheet instability, and include parameterizations for cliff failure and hydrofracturing of buttressing ice shelves (de Boer et al., 2015;Pollard et al., 2015;Gasson et al., 2016). However, such models require groundtruthing from the geological record.
The Pliocene (5.3 to 2.6 Ma) is a particularly relevant time period to study in light of future environmental change. Previous studies have suggested substantial changes in the Pliocene cryosphere, with largely ice-free conditions in the Northern Hemisphere (Haywood et al., 2016 and reference therein) and collapse of the West Antarctic Ice Sheet (WAIS) Pollard and DeConto, 2009). Estimates of global mean sea level (GMSL) of ∼22 ± 10 m higher than present during the warm Pliocene (Miller et al., 2012) may require not only collapse of the vulnerable polar ice in Greenland and West Antarctica, but also significant contributions from East Antarctica (see also Dutton et al., 2015 for discussion).
In the following, we present the first orbitally-resolved sediment provenance records from offshore of the Wilkes Subglacial Basin during the middle and late Pliocene. Our new data yield intriguing insights into the timescales and mechanisms of the equi-librium response of the EAIS to warmer than present environmental conditions in the geological past, and thereby provide important constraints for ice sheet modeling of future ice sheet behavior.

Age model at U1361A
Linear interpolation between paleomagnetic tie points is used to date U1361A Pliocene material following Tauxe et al. (2012) and Patterson et al. (2014) (Supplementary Table 1). The section studied is mostly continuous with one core gap between ∼3.6 to ∼3.33 Ma, which is thought to be related to the "super-glacial" M2 (Tauxe et al., 2012). We sampled for sediment provenance and biogenic silica concentrations (wt% BSi) at sub-orbital resolution (∼10 cm sampling corresponding to ∼2-4 kyr resolution).
Areas of high core disturbance from either drilling or bioturbation were avoided. Two of the high-resolution intervals studied, the Plio-Pleistocene boundary (47.55-50.35 mbsf, ∼2.50 Ma), and the  (Lisiecki and Raymo, 2005), with paleomagnetic tie points in U1361A (Tauxe et al., 2012) marked by grey dotted lines; (C) lithostratigraphy at U1361A (Escutia et al., 2011) with boxes marked around the intervals studied; (D, E) Nd and Sr isotopic compositions of the fine fraction of detrital sediment; representative error associated with Nd analyses is shown in black, while error associated with Sr measurements is smaller than the data symbols; (F) biogenic silica content (this study, Taylor-Silva and Riesselman, 2018); (G) Ba/Al XRF scan data Patterson et al., 2014); (H) K/Ti XRF scan data; (I) iceberg rafted debris mass accumulation rates (IBRD MAR; Patterson et al., 2014). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) late Pliocene (from 64.05-67.85 mbsf, ∼3.08 Ma), are constrained by at least one direct paleomagnetic datum within the depth interval analyzed for provenance, and hence are well constrained in their ages. Age uncertainty increases away from these age tie points due to the assumption of linear sedimentation rates, which does not capture changes over glacial-interglacial transitions. Bio-turbation within the core ranges from millimeters to centimeters (Escutia et al., 2011), which is less than our sampling resolution. Furthermore, all data presented were collected from the same core material, making our interpretations on timing and sequences of deglacial events in the Wilkes Subglacial Basin independent of the age model applied.

Sample categorization
Samples have been categorized into three groups based on previously published data: warmer conditions (i.e. interglacial intervals), colder conditions (i.e. glacial intervals), and transitional samples. Firstly, lithofacies analysis was used to classify the samples into two groups of diatom-rich/bearing mudstones or massive and laminated mudstones (Patterson, 2014). Samples from anomalous lithologies (e.g., rare sand laminae associated with mass flow deposition) were excluded from this classification. Secondly, average ratios of X-ray fluorescence (XRF) scan records (Ba/Al; Cook et al., 2013;Patterson et al., 2014) and results from shipboard physical property measurements (Gamma Ray Attenuation (GRA) Bulk Density, and Natural Gamma Radiation (NGR); Escutia et al., 2011) were used to define criteria for warmer and colder conditions. Sample categorization was determined based on the average division of these three variables (warmer conditions: Ba/Al > 0.231, GRA Bulk Density < 1.575, NGR < 32.3; colder conditions: Ba/Al < 0.171, GRA Bulk Density > 1.597, NGR > 34.9).
Where these conditions were not met, samples were classified as transitional. All chosen values are within the parameters constrained by previous studies at the Wilkes Subglacial Basin margin including high Ba/Al during warmer periods Patterson et al., 2014) and higher NGR during colder periods (NGR > ∼34; Tauxe et al., 2015).

Strontium and neodymium isotope analyses
For high resolution provenance analyses of the detrital fraction, residual sample splits from iceberg rafted debris mass accumulation rates (IBRD MAR) studies (>150 μm; Patterson et al., 2014) were utilized and combined with new samples wet sieved at <63 μm or <150 μm. In total 52 samples were selected for strontium (Sr) isotope analyses and 72 samples for neodymium (Nd) isotope analysis. A representative ∼0.5 g of sediment was taken from the homogenized fine-grained material and subjected to a sequential leaching procedure to remove biogenic carbonate using buffered acetic acid (see Biscaye (1965) for more details) and ferromanganese oxides and oxyhydroxides using a reductive solution of 0.02 M hydroxylamine hydrochloride (NH 2 OH) (Chester and Hughes, 1967;Cook et al., 2013Cook et al., , 2017. As in previous studies, no removal of biogenic opal was conducted (e.g. Cook et al., 2013). Approximately 50 mg of leached detrital sediment was digested on a hotplate using 2 ml of 27 M HF, 1 ml of 15 M HNO 3 and 0.8 ml of 20 M HClO 4 , followed by a secondary digestion step using 3 ml 15 M HNO 3 , before fluorides were removed with 3 ml 6 M HCl.
Ion chromatography was used to extract the target analytes (Sr and Nd) from the sample matrix. Samples were initially processed through a series of columns consisting of Tru-Spec resin (100-120 μm bead size), Ln-spec resin (50-100 μm bead size, 800 μl reservoir size) and Eichrom Sr-Spec resin (100-120 μm bead size). Later processing of samples used Biorad cation exchange resin (200-400 μm mesh), Ln-spec resin (50-100 μm bead size, 500 μl reservoir size) and Eichrom Sr-Spec resin (100-120 μm bead size). The change in procedure for separating REEs from the sample matrix in the first column resulted in improved Nd yields by preventing leakage of organics from Tru-Spec resin, subsequently affecting yields on Ln-Spec columns (see Struve et al., 2017 andLambelet et al., 2016 for details).
Neodymium isotope analyses were performed on a Nu Instruments multi collector inductively coupled plasma mass spectrometer (MC-ICP-MS) in the MAGIC laboratories at Imperial College London. Instrumental mass bias was corrected using the exponential law and a 146 Nd/ 144 Nd ratio of 0.7219. A correction for direct 144 Sm interference was applied, but all samples were significantly below the threshold determined for accurate correction (<0.1% of the 144 Nd signal). All reported sample 143 Nd/ 144 Nd ratios are corrected using bracketing standards and normalization to the recommended JNdi-1 value of 0.512115 (Tanaka et al., 2000). Measured average JNdi-1 values for individual analytical sessions are provided in Supplementary Table 2. Sixteen replicates of standard reference material BCR-2 were processed alongside unknown samples, and yielded an average 143 Nd/ 144 Nd ratio of 0.512639 ± 0.000022 (2 S.D., n = 46), identical to the recommended value of 0.512638 ± 0.000015 (2 S.D.; Weis et al., 2006). Procedural blanks were typically less than 85 pg, and always below 215 pg, which is negligible compared to the sample size processed (ranging between ∼200 and ∼1000 ng). Duplicate sample analysis (i.e., full procedural repeats from digestion through ion exchange chromatography to mass spectrometry) and re-analyzed samples (i.e. the same aliquot from ion exchange chromatography) always yielded results that agreed within error (see Supplementary Table 2).
Strontium isotopes were analyzed on a Thermo Scientific Triton thermal ionization mass spectrometer (TIMS) at the MAGIC laboratories at Imperial College London. Samples were loaded in 1 μl 6 M HCl onto degassed single tungsten filaments, followed by 1 μl of tantalum chloride activator. Instrumental mass bias was corrected using the exponential law and an 88 Sr/ 86 Sr ratio of 8.375, while interferences of 87 Rb were corrected using an 87 Rb/ 85 Rb ratio of 0.3860. Measurements were made in static mode using amplifier rotation. Repeated analysis of NBS987 standards over the duration of sample analysis yielded an average 87 Sr/ 86 Sr value of 0.710265 ± 0.000016 (2σ S.D.; n = 29). All reported 87 Sr/ 86 Sr ratios were corrected relative to the published value for NBS987 of 0.710252 ± 0.000013 (Weis et al., 2006) (Supplementary Table 2 Table 2), whereas total procedural duplicates of sediment showed a larger uncertainty (i.e. in the fourth decimal place). This variability is likely due to sample heterogeneity, and is comparable to what has been reported in other studies Tütken et al., 2002). Procedural blanks were typically less than 70 pg, representing considerably less than a 0.5% contribution to sample measurements.
In addition, tests were conducted on selected samples to ensure that using two different grain size fractions (i.e., <63 μm and <150 μm) of detritus yielded comparable results (see Supplementary Material).

Biogenic silica
Biogenic silica concentrations (wt% BSi) were measured on discrete samples using a molybdate blue spectrophotometric method modified from Strickland and Parsons (1972) and DeMaster (1981).
For each analysis ∼7 mg of dry, homogenized sediment was leached in 0.1 M NaOH at 85 • C, and sequential aliquots were collected after 2, 3, and 4 hours. After the final sampling hour, molybdate reagent and reducing solution were added to all samples, blanks, and standards to produce a colorimetric reaction, and absorbance of the 812 nm wavelength was measured using a Shimadzu UV-1800 spectrophotometer. Dissolved silica concentration of unknowns were calculated using a 10-point standard curve with known concentrations ranging from 0 μM to 1200 μM, and data from the three sampling hours were regressed following the method of DeMaster (1981) to calculate wt% BSi. Reproducibility was assessed by replication within each run and across runs. The average standard deviation of replicate measurements was 0.47%.

Provenance of Pliocene detrital sediments at IODP Site U1361A
Provenance tracing of detrital sediment is an effective tool for deciphering ice sheet histories (Licht and Hemming, 2017). The Wilkes Subglacial Basin is an ideal location for such studies due to the distinct geochemical signature of the lithological terranes in the hinterland Fig. 1b). Our new detrital sediment radiogenic isotope measurements provide evidence for gradual provenance changes across three distinct Pliocene glacialinterglacial cycles (Fig. 2), with consistent co-variation of Sr and Nd isotope fingerprints. In line with previous interpretations from low resolution records at the same site , transitions from cooler to warmer periods are characterized by shifts to less radiogenic (lower) 87 Sr/ 86 Sr ratios and more radiogenic (higher) ε Nd values (Fig. 2). As such, two distinct isotopic clusters are identified for the Pliocene detrital material at IODP Site U1361A: (1) ε Nd ranging from −12.6 to −15.5 and 87 Sr/ 86 Sr from 0.729 to 0.740, and (2) ε Nd = −9.0 to −11.0 and 87 Sr/ 86 Sr = 0.715-0.723. These denote (1) erosion of Paleozoic granites ( 87 Sr/ 86 Sr = 0.712 to 0.753, ε Nd = −9.7 to −19.8), exposed around the proximal Ninnis Glaciers (i.e. coastal outcrops), and (2) influence of an additional erosion source of inland material within the Wilkes Subglacial Basin (i.e. Ferrar Large Igneous Province and associated Beacon Supergroup (FLIP: 87 Sr/ 86 Sr = 0.709 to 0.719, ε Nd = −3.5 to −6.9)) (Fig. 3). The FLIP and Beacon Supergroup lithologies are located mostly inland and constitute the majority of the sedimentary infill in the Central Basin of the Wilkes Subglacial Basin (Ferraccioli et al., 2009;Frederick et al., 2016;Studinger et al., 2004) (Fig. 1b). Since such lithologies are currently covered by ice, significant erosion from this source is not detected in modern core top samples . Therefore, the appearance of this provenance signature in diatom-rich Pliocene sediments is interpreted to represent significant retreat of the ice margin into the Central Wilkes Subglacial Basin .

Sequence of Pliocene deglaciation events around the Wilkes Subglacial Basin
Our high resolution provenance study allows us, for the first time, to place variations in the source of the detrital material in direct context with other environmental changes documented in the same core, making the relative timing of events robust and independent of age model uncertainty. Firstly, we compare our new data to marine productivity proxies; continuous Ba/Al ratios Patterson et al., 2014) from XRF scans, and biogenic silica data (this study, see also Taylor-Silva and Riesselman, 2018) (Fig. 2). The shift in the geochemical provenance signal coincides closely with changes in both productivity records. Periods of climatic warmth (i.e. interglacial intervals) are characterized by biogenic-rich facies, with higher Ba/Al ratios (>0.21) and elevated biogenic silica content (8.5-25.4 wt%) (Fig. 4). Colder periods (i.e. glacial intervals) are characterized by biogenic-poor facies, lower Ba/Al ratios (<0.17) and lower biogenic  Cook et al., 2013). Analyses from samples categorized within this study as warmer conditions (i.e. interglacial intervals), colder conditions (i.e. glacial intervals), and transitional times are indicated by the color and symbol selection. Grey triangles denote previously published Pliocene data , which also represent warmer conditions, colder conditions and transitional climate states. (For interpretation of the colors in the figure, the reader is referred to the web version of this article.) silica content (5.2-15.4 wt%). The new biogenic silica data supports the interpretation made that diatom-rich/bearing mudstones with high Ba/Al ratios are indicative of times with higher productivity and hence warmer climatic periods Patterson et al., 2014).
When compared to iceberg rafted debris mass accumulation rates (IBRD MAR; Patterson et al., 2014), the timing of rapid shifts in sediment provenance and increase in productivity coincides with peak IBRD accumulation (Fig. 2). Maxima in IBRD have been interpreted to reflect accelerated calving of marine terminating glaciers , which would lead to a dynamic retreat of the ice sheet and subsequent re-stabilization of the ice margin inland, resulting in the observed provenance shifts. Large-scale ice retreat is therefore more likely to cause the observed provenance shifts than a switch in sediment delivery between ice streams. Marine productivity in the local Southern Ocean would subsequently increase as a result of the accompanying decrease in sea ice extent (Armbrecht et al., 2018;Taylor-Silva and Riesselman, 2018) and associated changes in nutrient delivery, allowing marine productivity to increase, as documented in the Ba/Al and biogenic silica records (Fig. 2).

Ice retreat and Southern Ocean biogeochemistry
We observe a significant negative correlation between XRF Ba/Al and K/Ti ratios in Pliocene marine sediments ( Fig. 5; R 2 = 0.85). K/Ti ratios indicate the provenance of terrigenous material (Monien et al., 2012), providing a continuous and independent provenance tracer to Nd and Sr isotopes. The inferred provenance shifts at IODP Site U1361A during interglacial intervals introduce a basaltic component to the source material, which carries lower K/Ti ratios, has substantially elevated Fe concentrations and has a higher susceptibility to weathering. Since Fe is a limiting nutrient in the modern Southern Ocean, the supply of highly soluble  Higher Ba/Al is interpreted as increased marine productivity, while lower K/Ti is indicative of an increased basaltic influence in the source sediment, interpreted to be caused by a retreated ice margin within the Wilkes Subglacial Basin.
reactive Fe in the form of labile (oxyhydr)oxides from glaciated terrains (Raiswell et al., 2006) could have played a major role in providing bioavailable trace metals and nutrients, driving Pliocene glacial-interglacial productivity changes.
However, marine productivity did not increase immediately at the time of increased IBRD accumulation (Fig. 2), suggesting that ocean fertilization around icebergs (Duprat et al., 2016) was not the dominant process driving enhanced interglacial productivity. Instead, our data supports the idea that bioavailable Fe trapped within multi-year sea ice (Geibert et al., 2010) could have been released during interglacials. Alternatively, or in addition, increased meltwater fluxes from the recently deglaciated terrains and/or enhanced subglacial erosion may have been important for nutrient delivery to the ocean. Statham et al. (2008) pointed out the large potential for glacially-derived dissolved and colloidal Fe to impact Southern Ocean primary productivity, and Hodson et al. (2017) showed that glacial meltwaters in the Antarctic Peninsula, where the source rock predominantly comprises young basaltic lithologies, contain significant amounts of bioavailable Fe. Additionally, in the Labrador Sea, offshore of the Greenland Ice Sheet, Arrigo et al. (2017) showed that the nutrients supplied in such meltwaters travel vast distances (>500 km) offshore of the coast. Sedimentary facies analysis from the Ross Sea furthermore indicates that the EAIS margin had a similar glacial regime to that of the modern Antarctic Peninsula and East Greenland during the Late Pliocene and Early Pliocene, respectively , which is consistent with enhanced glacial meltwater delivery at these times and hence enhanced localized nutrient delivery.

Timescales of ice retreat and the Pliocene history at the Wilkes Subglacial Basin ice margin
Constraining the history of the Antarctic ice sheet is important for improving our ability to project future rates of ice retreat in a warming world. However, obtaining rate information from paleorecords is inherently challenging because of bioturbation, relatively low sedimentation rates, and variable sediment deposition along glaciated margins. Nevertheless, the well-defined chronology and near continuous sediment deposition at U1361A over the Pliocene (Tauxe et al., 2012) allows linear sedimentation rates (0.85 to 5.49 cm/kyr) between direct datums to be utilized to extract information on timescales (Supplementary Table 1). As such, a first data-based estimate for natural timescales of ice retreat under elevated Pliocene temperatures is determined to be in the order of ∼3 to 7 kyrs. Timescales of deglaciation were estimated based on the distance between samples with isotopic compositions of colder conditions and the first sample representative of warmer conditions. This was found to be 20 cm, 10 cm and 20 cm (∼7, 3 and 5 kyr respectively) for the Plio-Pleistocene Boundary, Late Pliocene and Middle Pliocene intervals, respectively. The high fidelity orbital signature of the IBRD signal identified by Patterson et al. (2014) gives us confidence that we are assessing single eccentricity-paced (100 kyr) cycles in the time slices centered on 3.1 and 2.5 Ma, and a single obliquity cycle at 3.9 Ma (Fig. 2). Hence, our geological data provides support that ice margin retreat occurred on timescales similar to recent glacial terminations and modeling results from the Wilkes Subglacial Basin (i.e. <10 kyr; Pollard et al., 2015;Golledge et al., 2017b).
Our new record furthermore reveals similar provenance changes during all Pliocene intervals, providing evidence for continued ice retreat as global climate transitioned from Antarctic-dominated eccentricity control to a bipolar obliquity-paced state in the late Pliocene. While previous studies have suggested changes in ice dynamics around East Antarctica associated with late Pliocene global cooling and Antarctic sea ice expansion (Cook et al., 2014;Escutia et al., 2009;McKay et al., 2012;Orejola et al., 2014;Taylor-Silva and Riesselman, 2018), our new data demonstrates that the Wilkes Subglacial Basin was a highly sensitive ice-oceanatmosphere system in all climate states, despite significant Antarctic cooling at 3.3 Ma (McKay et al., 2012; Riesselman and Dunbar, 2013;Taylor-Silva and Riesselman, 2018), and the onset of Northern Hemisphere Glaciation at 2.7 Ma (Lisiecki and Raymo, 2005;Raymo et al., 2006). Throughout the late Pliocene, this region remained susceptible to ice retreat, warranting future investigations on its sensitivity in the Late Pleistocene under temperatures similar to those expected for the end of the 21st century.

Conclusions
Our results constitute the first orbitally-resolved provenance studies on one of the most vulnerable areas of the East Antarctic Ice Sheet under warmer environmental conditions and elevated atmospheric carbon dioxide concentrations. We demonstrate dynamic behavior of the ice margin at the mouth of the Wilkes Subglacial Basin throughout the Pliocene epoch, with ice retreat on timescales of a few thousand years, substantiating findings of recent modeling studies. Comparison of sediment provenance, marine productivity, and iceberg calving rates from the same samples indicates that the first deglacial response to orbitally-forced atmospheric and/or oceanic warming was an intensification in iceberg calving events (i.e. response of individual ice streams). This change triggered a retreat of the ice margin inland, resulting in a shift in the provenance signature and an increase in ocean productivity. Elevated productivity coincident with ice sheet retreat indicates that both reduced sea ice cover and glacial meltwater production could have played a role in local fertilization of the Southern Ocean during past warm periods. The Wilkes Subglacial Basin did not permanently stabilize within the Pliocene epoch, suggesting it may be more sensitive to melting than previously thought.

Acknowledgements
Provenance analysis was supported by a Kristian Gerhard Jebsen PhD Scholarship and NERC UK IODP grants (NE/H025162/1 and NE/H014144/1). Biogenic silica data was supported by a Royal Society of New Zealand Marsden FastStart grant (#UOO-1315) and a University of Otago PhD Scholarship. Support for sedimentology analysis was provided by the Royal Society of New Zealand Rutherford Discovery Fellowship (RDF-13-VUW-003). XRF work was supported by the Ministry of Science and Innovation Grant CTM2014-60451-C2-1-P co-financed by the European Regional Development Fund (FEDER). Samples were provided by the Integrated Ocean Drilling Program. The authors acknowledge two anonymous reviewers whose thoughtful comments helped improve the manuscript.

Appendix A. Supplementary material
Supplementary material related to this article can be found online at https://doi .org /10 .1016 /j .epsl .2018 .04 .054.