Ice sheet model simulations

Here the results of high-resolution ice sheet modelling experiments that investigate the dynamic response of an LGM-configuration Antarctic ice sheet to ocean forcing (Golledge et al. Reference Golledge, Fogwill, Mackintosh and Buckley2012, Reference Golledge, Levy, McKay, Fogwill, White, Graham, Smith, Hillenbrand, Licht, Denton, Ackert, Maas and Hall2013) are presented. The parallel ice sheet model (PISM) is used, a 3-D, thermomechanical, continental ice sheet model, constrained by published geological data that define lateral and vertical extents of the expanded Antarctic ice sheets around the time of the LGM (Golledge et al. Reference Golledge, Fogwill, Mackintosh and Buckley2012). The model combines the shallow-ice and shallow-shelf approximation equations across the entire domain. Therefore, the model is able to capture the dynamic behaviour within grounded ice of Antarctic ice sheets and simulate the drawdown of interior ice by ice streams at high resolution (5 km). In this study, the model uses proxy-based interpretations of oceanic (Lisiecki et al. Reference Lisiecki and Raymo2005, Imbrie et al. 2006) and atmospheric (Petit et al. Reference Petit, Jouzel and Raynaud1999) changes during the last glacial cycle and employs boundary distributions from modified Bedmap topography (Le Brocq et al. Reference Le Brocq, Payne and Vieli2010), temperature and precipitation fields from gridded datasets (Comiso Reference Comiso2000, Van de Berg Reference Van de Berg, van den Broeke, Reijmer and van Meijgaard2006), and a spatially varying geothermal heat flux interpolation (Shapiro & Ritzwoller Reference Shapiro and Ritzwoller2004).

The ice sheet model computes ice thickness and temperature changes, isostatic depression of topography, migration of grounding lines and the growth of ice shelves. Interaction between modelled ice shelves and their surrounding ocean is accounted for using a mass balance determination based on heat flux across the ice-water boundary. Our perturbation experiments use isochronous changes to oceanic heat flux and sea level values. The ice sheet model simulations are based on ocean-perturbation experiments in which the oceanic heat flux and sea level are isochronously increased from 30% to 100% of glacial to interglacial transition values, and by 25 m and 50 m, with respect to LGM values, representative of Holocene values. The response of the ice sheet is considered in terms of changes in velocity and ice thickness which together yield mass flux and, importantly, changes in ice flow direction.

Geological and geochronological constraints

To provide a temporal context for the modelled ice stream response, new and recalibrated existing geochronological data are used to reconstruct the geometric and temporal changes in the ice streams feeding the western part of the Weddell Sea following the LGM. The terrestrial record constrains altitudinal changes of the former ice stream surface, whereas offshore marine records, from radiocarbon dating of glaciomarine sediments overlying the subglacial deposits, constrain the lateral extent of the ice sheet in the WSE. Both are required to reconstruct the 3-D changes in ice stream geometry and investigate palaeo ice volume changes.

Marine geochronological constraints

The limited available analyses of marine sediment cores (including radiocarbon) from the outer and inner continental shelf of the WSE are used to constrain the lateral extent of the ice sheet, and provide the timing of ice sheet grounding line retreat and establishment of open water conditions (Larter et al. Reference Larter, Graham, Hillenbrand, Smith and Gales2012, Stolldorf et al. Reference Stolldorf, Schenke and Anderson2012, Hillenbrand et al. Reference Hillenbrand, Bentley, Stolldorf, Hein, Kuhn, Graham, Fogwill, Kristoffersen, Smith, Anderson, Larter, Melles, Hodgson, Mulvaney and Sugden2014). Whilst the existing marine chronology is open to interpretation, including possible reworking and the potentially significant changes in Antarctic marine radiocarbon reservoir effect over time (Hillenbrand et al. Reference Hillenbrand, Bentley, Stolldorf, Hein, Kuhn, Graham, Fogwill, Kristoffersen, Smith, Anderson, Larter, Melles, Hodgson, Mulvaney and Sugden2014), the ages are internally coherent, suggesting that they provide reliable constraints on the gradual retreat of the grounding line across the Weddell Sea. An important constraint from the outer continental shelf records grounding line retreat, and suggests open water conditions were established by c. 18.1 ka (Hillenbrand et al. Reference Hillenbrand, Bentley, Stolldorf, Hein, Kuhn, Graham, Fogwill, Kristoffersen, Smith, Anderson, Larter, Melles, Hodgson, Mulvaney and Sugden2014). Two further reliable critical constraints exist close to the sills of the extensive Thiel and Rutford cross-shelf troughs in the Weddell Sea. These record retreat in the eastern Weddell Sea at the head of the Thiel Trough before c. 8.3 ka, and at the head of the Rutford Trough in the western Weddell Sea at c. 5.3 ka (Fig. 1) (Hillenbrand et al. Reference Hillenbrand, Bentley, Stolldorf, Hein, Kuhn, Graham, Fogwill, Kristoffersen, Smith, Anderson, Larter, Melles, Hodgson, Mulvaney and Sugden2014).

Fig. 1 Weddell Sea embayment (WSE) indicating the sampling locations next to the Rutford and Institute ice streams. Ice sheet surface velocity data (Rignot et al. Reference Rignot, Mouginot and Scheuchl2011) highlight the locations of the major ice streams in light colours, and ice rises and slow moving regions in the WSE in darker blue. The sites of marine cores and associated minimum ages for grounding line retreat based upon marine radiocarbon ages (Hillenbrand et al. Reference Hillenbrand, Bentley, Stolldorf, Hein, Kuhn, Graham, Fogwill, Kristoffersen, Smith, Anderson, Larter, Melles, Hodgson, Mulvaney and Sugden2014) are also shown. F=Flower Hills, U=Union Glacier, P/M=Patriot and Marble hills.

Terrestrial geochronological constraints

Although the available marine radiocarbon data is sparse, a more comprehensive terrestrial record of ice stream surface changes is recorded on exposed mountains in the catchments of the Rutford and Institute ice streams. This study combines new ¹⁰Be and in situ ¹⁴C data that record changes of the Rutford Ice Stream with published in situ ¹⁰Be and ²⁶Al cosmogenic isotope data from the catchment of the Institute Ice Stream (Bentley et al. Reference Bentley, Fogwill, Le Brocq, Hubbard, Sugden, Dunai and Freeman2010, Fogwill et al. Reference Fogwill, Hein, Bentley and Sugden2012). Terrestrial ice surface elevations through time are constructed by measuring cosmogenic nuclides in erratics glacially transported from sites located in the catchments of the Rutford Ice Stream and Institute Ice Stream, as suggested by our LGM ice sheet flow model (Fig. 2). Glacial erratics sampled from steep exposed bedrock surfaces, in the Flower Hills, Union Glacier, and the Patriot and Marble hills (Bentley et al. Reference Bentley, Fogwill, Le Brocq, Hubbard, Sugden, Dunai and Freeman2010, Fogwill et al. Reference Fogwill, Hein, Bentley and Sugden2012) (Fig. 1), serve as ‘dipsticks’ that allow us to reconstruct past surface elevation changes in the catchments of the ice streams since the LGM.

Fig. 2 Simulated regional ice flux (upper panels), together with ice flow direction (white arrows) and ice sheet surface elevation of the Rutford and Institute ice streams (lower panel). a. Post-LGM conditions. b. Initial response to imposed ocean forcing leads to widespread acceleration of ice flow at principal outlets at c. 15 000 model years. c. Continued ice recession then leading to capture of the Institute Ice Stream by the Thiel Trough outlet during the late to mid Holocene. Ice flow vectors in the area of interest illustrate the change in flow direction taking place between time slices and red squares show the sample locations. F = Flower Hills, UG = Union Glacier, P/M = Patriot and Marble hills.

Samples were reduced to pure quartz at the University of Edinburgh cosmogenic nuclide laboratory and Lawrence Livermore National Laboratories Center for Accelerator Mass Spectrometry (LLNL-CAMS) following standard procedures (Kohl & Nishiizumi Reference Kohl and Nishiizumi1992, Ivy-Ochs Reference Ivy-Ochs1996, Stone Reference Stone2004). The ¹⁰Be ratios were measured by the AMS facility at LLNL and the Scottish Universities Environmental Research Centre (SUERC) (Xu et al. Reference Xu, Dougans, Freeman, Schnabel and Wilcken2010). Measurements were standardized to the NIST SRM-4325 Be standard material with a revised nominal ¹⁰Be/⁹Be ratio of 2.79x10^-11 (Nishiizumi et al. Reference Nishiizumi, Imamura, Caffee, Southon, Finkel and McAninch2007). Samples were corrected for the number of ¹⁰Be atoms in their associated blanks. Blanks were spiked with 250 μg ⁹Be carrier (Edinburgh) and 474 μg ⁹Be (LLNL-CAMS). The corresponding combined process and carrier blanks ¹⁰Be/⁹Be ratios range between 1.6–5.47x10^-15. Sample and blank ¹⁰Be/⁹Be analytical uncertainties and a 2.5% carrier addition uncertainty are propagated into the 1σ analytical uncertainty for nuclide concentrations.

A version of the CRONUS-Earth online age calculator was used to determine the ¹⁰Be exposure ages (Balco et al. Reference Balco, Stone, Lifton and Dunai2008), implementing the New Zealand ¹⁰Be production rate calibration dataset (Putnam et al. Reference Putnam, Schaefer, Barrell, Vandergoes, Denton, Kaplan, Finkel, Schwartz, Goehring and Kelley2010), that uses the recently revised ¹⁰Be half-life (1.387 Ma) (Chmeleff et al. Reference Chmeleff, von Blanckenburg, Kossert and Jakob2010, Korschinek et al. Reference Korschinek, Bergmaier and Faestermann2010), and Be isotope ratio standardization of Nishiizumi (Nishiizumi et al. Reference Nishiizumi, Imamura, Caffee, Southon, Finkel and McAninch2007). The use of this revised production rate and half-life change impact the apparent exposure ages, causing them to increase by c. 12% from those previously published (Bentley et al. Reference Bentley, Fogwill, Le Brocq, Hubbard, Sugden, Dunai and Freeman2010, Fogwill et al. Reference Fogwill, Hein, Bentley and Sugden2012) (Table I). Choice of production rate model and scaling is often a pragmatic one and is an ongoing subject of debate. Here the New Zealand calibration dataset was used to allow comparison with other recent Antarctic studies and in the absence of an Antarctic production rate calibration site. Exposure ages are reported based on the Lal/Stone scaling model for Antarctica; using the same calibration dataset, ages differ by 2–4% depending on the choice of scaling model (Balco et al. Reference Balco, Stone, Lifton and Dunai2008). The calculator uses sample thickness and density to standardize nuclide concentrations to the rock surface. The whole rock density is assumed to be 2.65–2.7 g cm^-3. No correction for periodic snow cover or for rock-surface erosion was included, as both of which are assumed to be negligible in these sites. An erosion rate of 0.0002 cm yr^-1 increases ages by c. 2%.

Table I ¹⁰Be cosmogenic isotope data from the Patriot and Marble hills recording changes in the Institute Ice Stream, and data from the Flower Hills and Union Glacier recording changes in the Rutford Ice Stream (Fogwill et al. Reference Fogwill, Hein, Bentley and Sugden2012).

^aRatio of the production rate at the shielded site to that for a 2π surface at the same location calculated using the CRONUS-Earth geometric shielding calculator version 1.1.

^bCalculated using 07KNSTD ¹⁰Be measurement standard and calibration with a reported ¹⁰Be/⁹Be ratio 2.85x10^-1240 or to the NIST standard with an assumed isotope ratio of 2.79x10^-11 and ¹⁰Be half-life 1.36 Ma (Chmeleff et al. Reference Chmeleff, von Blanckenburg, Kossert and Jakob2010, Korschinek et al. Reference Korschinek, Bergmaier and Faestermann2010).

^cModel exposure age assuming no inheritance, zero erosion, density 2.65–2.7 g cm^-3 and standard atmosphere calculated using the CRONUS-Earth ¹⁰Be-²⁶Al exposure age calculator (Balco et al. Reference Balco, Stone, Lifton and Dunai2008) version 2.2 using a constant production rate model and scaling scheme for spallation of Lal (1991)/Stone (2000). Ages based upon global production rate (P_global) and New Zealand production rate (P_NZ) accordingly.

^dEd=University of Edinburgh, SUERC=Scottish Universities Environmental Research Centre, CAMS-LLNL=Centre for Accelerator Mass Spectrometry-Lawrence Livermore National Laboratory.

Uniquely, this study also takes advantage of recent technological developments in the extraction and measurement of in situ radiocarbon (¹⁴C) from quartz (Hippe et al. Reference Hippe, Kober, Baur, Ruff, Wacker and Wieler2009, Reference Hippe, Kober, Wacker, Fahrni, Ivy-Ochs, Akçar, Schluchter and Wieler2013), a cosmogenic nuclide with a considerably shorter half-life than that of ¹⁰Be (¹⁰Be=1.36x10³ kyr, ¹⁴C=5.73 kyr). As the relatively short half-life of ¹⁴C means that in situ ¹⁴C acquired on exposure during interglacials decays if the sample is covered by ice during a subsequent glacial, the apparent ¹⁴C age reflects the true minimum exposure age of the sample. Crucially, the disparity between the ¹⁰Be and ¹⁴C data allows the potential influence of prior exposure or recycling in this setting to be assessed (Lifton et al. Reference Lifton, Jull and Quade2001, White et al. Reference White, Fülöp, Bishop, Mackintosh and Cook2011).

In situ ¹⁴C extraction was performed at ETH Zürich following a modified protocol (Hippe et al. Reference Hippe, Kober, Baur, Ruff, Wacker and Wieler2009, Reference Hippe, Kober, Wacker, Fahrni, Ivy-Ochs, Akçar, Schluchter and Wieler2013). Quartz aliquots of c. 5 g were preheated at c. 700°C to remove atmospheric ¹⁴C contamination followed by the extraction of in situ ¹⁴C during heating to 1550–1600°C for 2×2 hours. The collected CO₂ gas was split into two samples before AMS measurement due to large gas amounts. Samples were then measured with the MICADAS AMS system using the gas ion source (Ruff et al. Reference Ruff, Wacker, Gaggeler, Suter, Synal and Szidat2007, Synal et al. Reference Synal, Stocker and Suter2007, Wacker et al. Reference Wacker, Bonani, Friedrich, Hajdas, Kromer, Nemec, Ruff, Suter, Synal and Vockenhuber2010). The number of ¹⁴C atoms obtained for both splits were summed prior to subtraction of the long-term average processing blank of (3.15±1.19)×10^{4 14}C atoms (± 1 standard deviation, n=24).

¹⁴C/¹⁰Be multi-isotope analysis

For this study, in situ ¹⁴C exposure ages were calculated with a sea level, high latitude (SLHL) spallogenic production rate of 11.40±0.9 at g^-1 y^-1 (Schimmelpfennig et al. Reference Schimmelpfennig, Schaefer, Akçar, Koffman, Ivy-Ochs, Schwartz, Finkel, Zimmerman and Schlüchter2014). As with ¹⁰Be, the production rate was scaled to altitude and latitude according to the scaling scheme of Lal/Stone. The contribution due to muon production was calculated using the freely accessible MATLAB code of the CRONUS-Earth online calculator (http.//hess.ess.washington.edu/math/al _be_v2/P_mu_total) (Balco et al. Reference Balco, Stone, Lifton and Dunai2008). In order to allow muon scaling for in situ ¹⁴C, parameters were adjusted based on the cross sections for ¹⁴C (Heisinger et al. Reference Heisinger, Lal, Jull, Kubik, Ivy-Ochs, Knie and Nolte2002a, Reference Heisinger, Lal, Jull, Kubik, Ivy-Ochs, Neumaier, Knie, Lazarev and Nolte2002b), and corrections for sample thickness and topographical shielding were applied on spallogenic production only.

Combined ¹⁴C and ¹⁰Be analysis is applied to sites in the Flower Hills and Union Glacier (in the catchment of the Rutford Ice Stream) which display a high percentage of anomalously ‘old’ apparent ¹⁰Be exposure ages (Table I). The disparity between the in situ ¹⁴C and ¹⁰Be data demonstrates that the samples have experienced a complicated exposure history, suggesting either that the cosmogenic nuclide inventories of the erratics were not fully reset by glacial erosion prior to deposition, or that following initial deposition they underwent periods of exposure at different altitudes and/or cover by cold-based ice (White et al. Reference White, Fülöp, Bishop, Mackintosh and Cook2011). Using the measured concentrations of both ¹⁰Be and ¹⁴C, an iterative model was constructed to calculate the maximum and minimum periods of ice cover each sample could have undergone to explain the differing nuclide concentrations. These periods were then compared with the equivalent periods of ice cover implied by the eustatic sea level data, following a similar approach to studies of the Fennoscandian Ice Sheet (Fabel et al. Reference Fabel, Stroeven, Harbor, Kleman, Elmore and Fink2002).

Whilst only a one-way test, this assumes that the erratics experienced periods of ice cover subsequent to initial deposition at any point when sea level was lower (and ice volume greater) than it was at the point of re-exposure given by the ¹⁴C apparent exposure age (Table II). Using the measured, minimum and maximum ¹⁴C exposure ages of each sample (given by the external errors within the process, to allow for comparability with the independently dated sea level curve), three scenarios for each sample were created, under which all samples apart from UG-27 are shown to have experienced one period of extended ice cover following initial deposition, followed by subsequent re-exposure (Fig. 3). Sample UG-27 has a complex nuclide inventory, possibly reconcilable by either a single or multiple pre-exposure event at a higher altitude than present. For the remaining samples, a scenario is identified that agrees with the periods of ice cover stipulated by the cosmogenic nuclide data and the sea level reconstruction: a single extended period of ice cover following initial deposition, and subsequent re-exposure during the last deglaciation (Fig. 3). This suggests that the disparity between the ¹⁴C/¹⁰Be is probably a result of cover by cold-based ice, and demonstrates that the use of another isotope paired alongside ¹⁴C can provide insights into the depositional history of the samples, allowing for a more confident interpretation of the surface trajectory of the Rutford Ice Stream.

Fig. 3 Modelled relationship between ¹⁰Be/¹⁴C isotope concentrations, time and sea level used as a proxy for global ice volume (Imbrie & McIntyre Reference Imbrie and McIntyre2006) for samples FLO/18/CJF and UG16. Proposed periods of sample exposure are defined by the grey boxes. The altitude and apparent exposure ages based upon the measured ¹⁰Be and ¹⁴C inventories of the samples are noted. The inset photo shows sample FLO/18/CJF, a quartzite erratic on striated agrilite bedrock typical of the samples analysed.

Table II ¹⁴C cosmogenic isotope data from the Flower Hills and Union Glacier.

^aNormalized to δ¹³C of -25‰_VPDB and AD 1950.

^bCalculated after eq. (1) in Hippe et al. Reference Hippe, Kober, Baur, Ruff, Wacker and Wieler2009.

^cBlank corrected; calculated after eq. (2) in Hippe et al. Reference Hippe, Kober, Baur, Ruff, Wacker and Wieler2009.

^dCorrected for sample thickness and topographical shielding; see Table I for correction factors.

^e Exposure age assuming no inheritance, zero erosion, density 2.7 g cm^-3, with a constant production rate and scaling scheme for spallation of Lal (1991)/Stone (2000).