An 83 000-year-old ice core from Roosevelt Island, Ross Sea, Antarctica

In 2013 an ice core was recovered from Roosevelt Island, an ice dome between two submarine troughs carved by paleo-ice-streams in the Ross Sea, Antarctica. The ice core is part of the Roosevelt Island Climate Evolution (RICE) project and provides new information about the past configuration of the West Antarctic Ice Sheet (WAIS) and its retreat during the last deglaciation. In this work we present the RICE17 chronology, which establishes the depth–age relationship for the top 754 m of the 763 m core. RICE17 is a composite chronology combining annual layer interpretations for 0–343 m (Winstrup et al., 2019) with new estimates for gas and ice ages based on synchronization of CH₄ and δ¹⁸O_atm records to corresponding records from the WAIS Divide ice core and by modeling of the gas age–ice age difference.

Novel aspects of this work include the following: (1) an automated algorithm for multiproxy stratigraphic synchronization of high-resolution gas records; (2) synchronization using centennial-scale variations in methane for pre-anthropogenic time periods (60–720 m, 1971 CE to 30 ka), a strategy applicable for future ice cores; and (3) the observation of a continuous climate record back to ∼65 ka providing evidence that the Roosevelt Island Ice Dome was a constant feature throughout the last glacial period.

1 Introduction

The Roosevelt Island Climate Evolution (RICE) project seeks to combine geophysical measurements with climate information from a well-dated ice core to improve constraints of the glacial history of the eastern Ross Sea (Conway et al., 1999). With this motivation, the RICE project drilled and recovered a 763 m long ice core from Roosevelt Island. Here, we first present new data sets from the RICE ice core and then describe the development of the RICE17 age-scale, a composite age-scale which combines annual layer counts (Winstrup et al., 2019) with new age constraints and spans the top 753.75 m.

The most precise chronologies for ice cores have been constructed by combining absolute age markers, for example volcanic ash layers, with annual layers counts from visual stratigraphy and from variations in chemical concentrations, stable isotope ratios of the ice, and electrical properties. Examples include the Greenlandic ice cores composite chronology (GICC05) which extends to 60 ka (Svensson et al., 2008) and the Antarctic WAIS Divide WD2014 chronology through the last 31 ka (Sigl et al., 2016). For the RICE ice core this strategy was only applicable above 165.02 m, the depth of the oldest absolute age marker, at 1251 CE (Winstrup et al., 2019). Annual layer interpretations were extended to 343.7 m (2649 years BP; all ages are reported as years before present, BP, where present is defined as 1950 CE; Winstrup et al., 2019). Below this depth, annual layers became too thin to reliably interpret the seasonal signals.

For the 343.7–754 m section of the RICE ice core, the best age constraints are from measurements of methane and the isotopic composition of molecular oxygen (δ¹⁸O_atm). Variations in these parameters are globally synchronous (Bender et al., 1994; Blunier et al., 1998), meaning that, to a first order, different ice cores will record the same atmospheric history. We establish the gas-phase depth–age relationship by matching the new RICE records to corresponding records from other ice cores with established chronologies. From the gas age-scale, an ice-phase age-scale was obtained by modeling the ice age–gas age offset (Δage). Δage is constrained by measurements of the isotopic composition of molecular nitrogen (δ¹⁵N-N₂) and estimates of temperature based on δD of ice. For overlapping depths, the gas-based ice age estimates are in excellent agreement with the annual-layer-counted age estimates. A final age-scale for the entire core, named RICE17, is presented and is a composite of the annual-layer-based age-scale for 0–343.7 m depth and the gas-based age-scale for deeper sections of the core.

The approach used for RICE17 is not unique but includes several refinements that improve the dating accuracy. Primarily, the RICE17 chronology benefited from the availability of high-resolution gas records from both the RICE ice core and the WAIS Divide ice core. Centennial-scale variability in methane was well captured for the last 30 ka in both data sets. Additionally, an automated matching routine, adapted from Huybers and Wunsch (2004) and novel in the application to ice cores, simultaneously matched methane and δ¹⁸O_atm records for ice up to 30 ka (0–720 m). The multiproxy approach to synchronization leads to strengthened age estimates. Using an automated routine resulted in a more objective match of the RICE gas records to the reference records and a Monte Carlo-based estimate of age uncertainty.

RICE17 is a continuous chronology with age monotonically increasing with depth until at least 64.6 ka. At 746.00 m, folding of the ice is observed, and the age-scale is segmented and no longer continuous. The oldest ice near the bottom of the core has a minimum age of 83 ka. We use these observations, with support from measurements of total air content (TAC), to interpret that the Roosevelt Island ice dome was stable throughout the last glacial maximum.

2 The Roosevelt Island ice core and glaciological history

Understanding the stability of the West Antarctic Ice Sheet (WAIS) is important for predicting future sea level rise (Church et al., 2013; Jevrejeva et al., 2014; Golledge et al., 2014; Pollard et al., 2015; DeConto and Pollard, 2016). Much of the ice sheet is grounded below sea level with the bed deepening towards the center of the ice sheet. This configuration is thought to be unstable and prone to rapid disintegration due to physical forces related to buoyancy (Hughes, 1973; Weertman, 1974; Schoof, 2007; Feldmann and Levermann, 2015), vulnerability to undercutting by “warm” subsurface currents (Robin and Adie, 1964; Shepherd et al., 2004), and ice cliff instability (DeConto and Pollard, 2016). Vulnerability of WAIS to future warming can be assessed by investigating how it has responded to different climate regimes in the past. Unfortunately, geologic evidence of the past size and extent of Antarctic ice sheets is spatially sparse, tends to have large chronological uncertainty, and is sometimes contradictory (Whitehouse et al., 2012; Anderson et al., 2014; Bentley et al., 2014; Clark and Tarasov, 2014; Halberstadt et al., 2016; McKay et al., 2016).

The Ross Embayment is the largest drainage of WAIS, both in terms of area and mass loss (Halberstadt et al., 2016). Historically, two scenarios have been proposed for the configuration of WAIS in the Ross Embayment during the Last Glacial Maximum (LGM) (Stuiver et al., 1981). In the “maximum scenario,” a thick and grounded ice sheet in the Ross Sea extended to the continental shelf break (Denton et al., 1989). Details of this scenario are supported by geomorphic features including grounding-zone wedges, which form at the terminus of marine-based glaciers (Shipp et al., 1999), and over-compressed diamictons, which are the result of thick overlying ice (Anderson et al., 1984, 1992). Evidence of high stands in the Transantarctic Mountains and the islands of the western Ross Sea (Denton and Marchant, 2000) as well as cosmogenic exposure dates on nunataks in Marie Byrd Land (Stone et al., 2003) also support this idea. In an alternate scenario, Denton et al. (1989) proposed that grounded ice in the Ross Sea was kept thin by fast-flowing ice streams. In this scenario, the retreat of WAIS during the last deglaciation may not have contributed significantly to sea level change. Studies of ice cores from Byrd Station (Steig et al., 2001) and Siple Dome (Waddington et al., 2005), glacial modeling (Parizek and Alley, 2004), and cosmogenic exposure dates from the Ohio Range (Ackert et al., 1999, 2007) all support this “minimal scenario.”

Roosevelt Island is an ice dome located in the eastern Ross Sea in West Antarctica (79.36^∘ S, 161.71^∘ W; elev. 550 m above sea level, Fig. 1). It is grounded on a submarine plateau (∼200 m b.s.l.) dividing the Whales Deep and Little America basins (Fig. 1). During the LGM these troughs were presumably occupied by the extension of the modern Bindschadler and MacAyeal ice streams (Ice Streams D and E, respectively) (Shipp et al., 1999), and Roosevelt Island would have been located along the main ice flow of WAIS.

Figure 1 Map of bedrock elevation in the Ross Sea Sector of Antarctica (referenced to WGS84 datum) (Fretwell et al., 2013). Gray dashed lines indicate ice sheet grounding lines and ice margins. Locations of the RICE (Roosevelt Island), Siple Dome, Byrd Station, WAIS Divide, and Taylor Dome ice cores are marked with black stars.

The 763 m long RICE ice core was drilled to bedrock near the summit of Roosevelt Island. In addition to the main deep core, a shallow core was drilled to 20 m, and several snow pits were sampled to understand recent climate (Bertler et al., 2018). Local mean annual air temperature is −23.5 ^∘C, and annual snow accumulation is estimated to be ∼0.22 m ice equivalent, based on annual layers in snow pits (Bertler et al., 2018). A cooler estimate of modern temperature, $-\mathrm{27.4}\pm \mathrm{2.4}$ ^∘C, based on ERA interim data from 1979 to 2012 was presented in Bertler et al. (2018), but this estimate is cooler than borehole thermometry measurements and previously published estimates for Roosevelt Island (Herron and Langway, 1980; Conway et al., 1999; Martín et al., 2006) and does not provide a good fit to the density profile in the firn model. Other estimates of recent accumulation at the RICE drill site range from 0.18 to 0.27 m ice per year, depending on method and time period (Winstrup et al., 2019; Bertler et al., 2018; Herron and Langway, 1980; Conway et al., 1999; Kingslake et al., 2014).

3 New data sets from the RICE ice core

3.1 Methane concentrations

The RICE discrete methane record was measured at Oregon State University (OSU) following methods described by Mitchell et al. (2011, 2013) with updates described in Appendix A. A total of 702 samples were measured at 583 distinct depths between 60 and 753 m (Fig. 2a, f). Samples from 406 depths were measured between 60 and 670 m, dating from ∼1970 CE to 11.87 ka, with a mean sample spacing of 28.75 years. Between 670 and 718.13 m the record spans 11.87 to 29.9 ka; 96 samples in this interval provide age resolution of 189 years. Age resolution decreases significantly for deeper depths. The interval from 718.53 to 746.00 m corresponds to 30.1–64.6 ka with a mean resolution of 548 years. The deepest dated ice is at 752.95 m with an age of 83 ka (±2 ka).

Figure 2 Data from the RICE ice core. Left panels (a–e) (60–665 m) cover the last 11.26 ka; right panels (f–j) (665–760 m) cover measurements from 11.26 to 83 ka and measurements from the 9 m of ice below the dated section. (a, f) Continuous methane measurements, gray, between 0 and 726 m depth are calibrated to discrete methane measurements, black. Beyond 726 m depth, raw CFA methane measurements are plotted. (b, g) Continuously measured δD smoothed with a 500-year LOESS filter and subsampled at 50-year spacing. (c, h) δ¹⁸O_atm measurements are corrected for gravitational enrichment in the firn layer using δ¹⁵N-N₂ (d, i). (e, j) Total air content (TAC) was measured in conjunction with discrete methane.

Methane was also measured continuously with a laser spectroscopy technique (Stowasser et al., 2012; Rhodes et al., 2013) during two separate continuous flow analysis (CFA) campaigns at GNS Science (Gracefield, New Zealand) in 2013 and 2014 (Pyne et al., 2018). The CFA methane record was affected by variability in airflow to the measurement instrument and fractures that allowed drill fluid and modern air into the melt head. Exclusion of these artifacts caused significant gaps in the record, particularly at depths below 676 m (12.6 ka). We consider the CFA record to be a supplement to the more robust but lower-resolution discrete data. Between 29.9 and 59.1 ka, the CFA record is critical for establishing age control. Below ∼746 m (64.6 ka), the CFA data are difficult to interpret because of gaps in the record, uncertainty in measurement depth, and uncertainty in methane calibration (Appendix B).

Several anomalously high discrete methane measurements appear between 44.6 and 50.9 ka (729.05–736.05 m) and below 64.6 ka (746 m) (Fig. 3a). In the former interval, methane is enriched by ∼30 ppb compared with the WAIS Divide CFA data (Rhodes et al., 2015). This is likely due to healed fractures that include modern air (Aydin et al., 2010). All RICE samples deeper than 500 m were visually inspected for fractures as they were prepared for measurement and for drill fluid during air extraction. In fractured ice it was common to see drill fluid in the flask as the samples were melted. However, neither observation was a strong indicator that a sample would have an elevated methane concentration. None of these high-concentration results were rejected.

Figure 3 Data from the RICE ice core plotted on the RICE17 age-scale. (a) RICE methane and (c) δ¹⁸O_atm records (gray) shown in comparison to target records from the WAIS Divide ice core on the WD2014 age-scale (Rhodes et al., 2015; Buizert et al., 2015b) (red) and NGRIP ice core (Baumgartner et al., 2014; Landais et al., 2007) on a modified GICC05modelext chronology (Wolff et al., 2010) (green). Solid triangles above panel (a) are gas age constraints from a Monte Carlo analysis; open diamonds are prior ACPs from visual matching (see text). (b) δD and (d) δ¹⁵N-N₂ used in a firn densification model to estimate the ice age–gas age offset. (e) RICE TAC measurements. (f) Mean annual layer thickness calculated from gas age control points adjusted for Δage.

3.2 Total air content

Total air content (TAC) is defined as the amount trapped air in a sample, reported here in units of cubic centimeters of air at STP per gram of ice (cm³ g⁻¹). TAC (Fig. 2e, j) was measured at OSU as part of the methane concentration measurement following methodology of Mitchell et al. (2015). TAC is influenced by accumulation and temperature, total summer insolation (Raynaud et al., 2007), thermal gradients in the firn from multiannual climate trends, and surface air pressure (Martinerie et al., 1992; Raynaud and Whillans, 1982).

Preparation of ice samples causes loss of air from bubbles which intersect the surface of the sample, an effect known as the cut-bubble effect (Martinerie et al., 1990). The amount of air lost due to the cut-bubble effect depends on the geometry of the sample and the bubble size. For bubble size, we assume that the bubble diameter at the bubble-close-off depth is the same as at the bubble-close-off depth at WAIS Divide and that the diameter shrinks exponentially to zero at the bubble–clathrate-transition depth (763 m). Ice samples were cut to uniform shapes to limit variability in TAC related to the cut-bubble effect so that TAC can be directly compared between samples. However, inclusion of some fractures was often unavoidable, and their contribution to gas loss is potentially large. TAC data were rejected when gas loss was believed to greatly impact the results, such as in samples with visible fractures or samples which consisted of multiple pieces. Of the 706 samples measured at OSU for TAC, 165 results were rejected based upon visual inspection of the sample. Many of these came from the 670 to 752.95 m (11.7–83 ka) interval where only 58 of 177 TAC measurements are considered reliable. Reproducibility of replicate TAC measurements is 0.7 % after application of a correction for the cut-bubble effect (reproducibility of 0.6 % without correction). Estimate of relative uncertainty in TAC measurements is 28 % with the largest contributions to the uncertainty from the shape of the sample and the assumed bubble diameter used for the cut-bubble correction.

The TAC record from the RICE ice core without the cut-bubble correction appears remarkably consistent (age weighted mean =0.1182 cm³ g⁻¹; age weighted standard deviation =0.0023 cm³ g⁻¹; n=410). After application of the cut-bubble correction, the TAC record shows a trend to higher values at more recent ages (Fig. 3e). This trend reflects the decrease in bubble size with depth and therefore a larger cut-bubble correction for shallower samples. Although this trend has a physical explanation, it is not statistically significant due to the uncertainty in the cut-bubble correction. We advise caution when interpreting the absolute values of TAC from the RICE ice core due to the potential artifacts caused by gas loss through fractures and the large uncertainty in the cut-bubble correction.

3.3 δ¹⁸O_atm and δ¹⁵N-N₂

δ¹⁸O_atm and δ¹⁵N-N₂ were measured on samples adjacent to the discrete methane samples (Fig. 2c, h and d, i, respectively). Analysis was conducted at Scripps Institution of Oceanography following Petrenko et al. (2006) and Severinghaus et al. (2009). Pooled standard deviation of replicate measurements is 0.006 ‰ for δ¹⁸O_atm and 0.0027 ‰ for δ¹⁵N-N₂ (both scales are relative to modern atmospheric composition).

Variations in δ¹⁸O_atm are primarily caused by changes in location and intensity of low-latitude rainfall that affect the δ¹⁸O of leaf water used in photosynthesis (Severinghaus et al., 2009; Landais et al., 2007, 2010; Seltzer et al., 2017) and changes in seawater δ¹⁸O caused by ice sheets on glacial cycles (Horibe et al., 1985; Bender et al., 1985, 1994; Sowers et al., 1993). Importantly, these variations are well known in independently dated ice cores, and the atmosphere is well mixed on the relevant timescales, so δ¹⁸O_atm variability may be used as a chronostratigraphic marker (Bender et al., 1994). These variations are well sampled in the RICE record until ∼64.6 ka (746.00 m), beyond which the chronology is no longer continuous.

Atmospheric N₂ is isotopically stable over very long timescales (Hattori, 1983; Sowers et al., 1989; Schwander, 1989). Variability in δ¹⁵N-N₂ in the ice core records primarily reflects changes in gravitational fractionation and thermal fractionation within the firn due to changes in surface temperature and accumulation rate (Schwander, 1989; Sowers et al., 1989). In Antarctica, temperature changes gradually, and the effect of thermal separation is minimal.

3.4 δD

The RICE δD record (Fig. 2b, g) was measured using off-axis integrated cavity output spectroscopy (OA-ICOS) using a Los Gatos Research (LGR) isotope water analyzer (IWA-35EP) in combination with a modified LGR water vapor isotope standard source (WVISS) evaporation unit (Emanuelsson et al., 2015). The data were calibrated using a time-weighted, two-point linear calibration with two standards. The cumulative uncertainty is calculated as the quadrature sum of Allan variance error, scatter over the averaging (15 s) interval (error of the variance) and calibration error (error of the mean). Final mean total uncertainty for the RICE δD record is 0.74 ‰ (Keller et al., 2018). The record shown in Fig. 3b has been smoothed with a 500-year low-pass filter (LOESS filter) and sampled at 50-year intervals. The smoothed δD record is provided in the Supplement to this paper.

4 Strategy for developing the RICE17 chronology

In this section, we describe how the RICE17 age-scale is developed. We start with age constraints in the gas-phase based on well-recorded variations in δ¹⁸O_atm and methane. We then discuss how the gas age-scale is used to estimate the gas-phase ice-phase age offset and thus the ice age-scale. Finally, we compare the ice age-scale to the age-scale derived from annual layer interpretations.

All age control points (ACPs), additional control points from the synchronization routine (floating control points, discussed below), and gas and ice ages interpolated to more closely spaced depths are provided in the Supplement.

4.1 An automated matching algorithm for synchronizing ice core records: 0–30 ka

An automated matching algorithm was adapted from Huybers and Wunsch (2004) to synchronize between gas records. Prior ACPs, shown as open triangles in Figs. 3 and 4, all correspond to well-defined increases or decreases in either methane or δ¹⁸O_atm. We assume that the uncertainty (2 standard deviations) for an ACP corresponds to the time elapsed between 25 % and 75 % of the change in either methane or δ¹⁸O_atm (Table 1). The goal of the routine is to iteratively adjust the interpolation between ACPs to improve a goodness of fit parameter by following steps 1–9 described below.

Figure 4 (a) RICE methane and (b) δ¹⁸O_atm records matched to those from the WAIS Divide ice core (Rhodes et al., 2015; Buizert et al., 2015b) (red lines) for the last 11.7 ka. FCPs from Monte Carlo routine are shown as solid triangles; prior constraints are shown as open diamonds. Lower panels (c, d) show the distribution of the gas ages for two particular depths, 459.05 m (4675 years BP; c) and 621.28 m (9012 years BP; d), resulting from the Monte Carlo analysis. The final age of these depths, resulting from the best realization, are shown as vertical gray lines.

Table 1 Prior gas age constraints and uncertainty are based on matching of features in the atmospheric history of methane and δ¹⁸O_atm. Name or description of feature are given in the notes and the primary parameter used to identify the feature is provided as the identifying variable. Age uncertainty is related to the duration of the feature; see text.

“Goodness of fit” is calculated as the χ² value comparing the normalized methane record from the RICE ice core to expectant values interpolated from the WAIS Divide ice core plus the analogous χ² value comparing the normalized δ¹⁸O_atm records. The algorithm can accept additional parameters for synchronization, such as CO₂ or N₂O, if available. In this analysis we normalized both the methane and δ¹⁸O_atm records to have the same mean (5) and variance (1) in order to equally weight their χ² values and their importance to the synchronization.

A realization starts by randomly perturbing the age of ACPs within their prescribed uncertainty (Table 1) to define an initial depth–age scale. The perturbed ACPs remain fixed throughout the realization. ACPs are given the subscript k, where k=1 is the shallowest (youngest) ACP. The records are then broken into N subsections which are distributed between ACPs. Subsections are designated with a subscript i, where $i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots }N$. Floating control points (FCPs) are defined as the bounding depths (ages) of these subsections, which will include the initial ACPs. Prior to any perturbations, the durations of the subsections (Δt) are approximately the same. We have chosen Δt to roughly match the recurrence interval of variations in methane in the reference record and so that each subsection contains about five methane and five δ¹⁸O_atm samples. The following steps are repeated to optimize goodness of fit:

1. Random scaling factors (p_i), which perturb the durations of the subsections, are drawn for each subsection from a normal distribution of μ=1 and σ₁=0.25.

   $$\begin{array}{}\text{(1)}& \mathrm{\Delta }{t}_i^{\prime }=\mathrm{\Delta }{t}_i\cdot p_i.\end{array}$$

2. Because random perturbations will change the length of time between the initial ACPs, we apply a second scaling so the perturbed chronology remains consistent with the prior ACPs. In this case, ACP_k and ACP_k+1 are the nearest gas ACPs bounding subsection i, respectively. ACP^′ is the perturbed age of the gas age constraint after step 1, and $\mathrm{\Delta }{t}_i^{\ast }$ is the duration of the subsection after the second scaling:

   $$\begin{array}{}\text{(2)}& \mathrm{\Delta }{t}_i^{\ast }=\mathrm{\Delta }{t}_i^{\prime }\cdot \left({\displaystyle \frac{{\text{ACP}}_{k+\mathrm{1}}-{\text{ACP}}_k}{{\text{ACP}}_{k+\mathrm{1}}^{\prime }-{\text{ACP}}_k^{\prime }}}\right).\end{array}$$

3. Mean “accumulation rate” of each subsection (${\stackrel{\mathrm{‾}}{A}}_i$) is calculated following Nye (1963), which assumes that the thickness of annual layers (${\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_i$) is the product of their original thickness and their relative depth ($\frac{{\stackrel{\mathrm{‾}}{z}}_i}{H}$):

   $$\begin{array}{}\text{(3)}& {\stackrel{\mathrm{‾}}{A}}_i={\stackrel{\mathrm{‾}}{\mathit{\lambda }}}_i\cdot {\left(\mathrm{1}-{\displaystyle \frac{{\stackrel{\mathrm{‾}}{z}}_i}{H}}\right)}^{-\mathrm{1}},\end{array}$$

   where ${\stackrel{\mathrm{‾}}{z}}_i$ is the mid-depth of the subsection i, H is the thickness of the ice sheet, and both ${\stackrel{\mathrm{‾}}{z}}_i$ and H are in ice equivalent units.

   Because the assumptions about thinning from Nye (1963) are too simple to describe the flow conditions at Roosevelt Island, we do not consider this result to be representative of the true accumulation history. The assumption is necessary for the purposes of interpolating age versus depth to ensure that annual layer thickness decreases smoothly between FCPs. An accumulation history that we believe is more accurate is calculated below with an alternative method, using a dynamic firn model (Sect. 4.4).

4. A perturbed chronology is only accepted if the following are true:

   • a.

     Intervals between FCPs ($\mathrm{\Delta }{t}_i^{\ast }$) are within a factor of 10 of the initial durations,

     $$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{10}}}\cdot {\displaystyle \frac{\left(t_{\mathrm{end}}-t_{\mathrm{1}}\right)}{N}}<\mathrm{\Delta }{t}_i^{\ast }<\mathrm{10}\cdot {\displaystyle \frac{\left(t_{\mathrm{end}}-t_{\mathrm{1}}\right)}{N}}.\end{array}$$
   • b.

     Mean accumulation rates inferred in step 3 are realistic (${\stackrel{\mathrm{‾}}{A}}_i>$0 cm ice equivalent per year; ${\stackrel{\mathrm{‾}}{A}}_{max}<\mathrm{75}$ cm ice equivalent per year), and

   • c.

     Mean accumulation rates of adjacent subsections are within a factor of 2 of each other.

   These conditions provide loose restrictions for continuity in the depth–age relationship although they may not be physically realistic for the site.

5. Ages for all RICE sample depths are interpolated from the perturbed FCPs (Step 3), assuming piecewise constant accumulation between FCPs and a linearly decreasing thinning function (Nye, 1963).

6. Goodness of fit is calculated on the perturbed age-scale (goodness of fit = ${\mathit{\chi }}_{{\mathrm{CH}}_{\mathrm{4}}}^{\mathrm{2}}+{\mathit{\chi }}_{{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{atm}}}^{\mathrm{2}}$).

7. When a perturbation improves the goodness of fit, that iteration becomes the base for subsequent perturbations.

8. The cycle is repeated until 20 sequential perturbations fail to improve the goodness of fit in step 6.

9. Starting with the FCPs from step 8, the above steps (1–8) are repeated 13 times, reducing the size of perturbations in step 1 from ${\mathit{\sigma }}_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{4}}$ to ${\mathit{\sigma }}_{\mathrm{13}}=\frac{\mathrm{1}}{\mathrm{128}}$.

Table 2 Model Parameters used for optimized correlation routine. Code for routine can be found in the Supplement.

A Monte Carlo analysis repeats these steps (Steps 1–9) 1000 times, initiated from a different prior depth–age relationship by randomly perturbing ACPs within their age uncertainty (Table 1). Parameters used for the synchronization are given in Table 2 for the 60.05–670.13 m and 670.13–719.30 m intervals. Code for the synchronization routine is provided in the Supplement.

We choose the best realization for gas age constraints for the RICE17 gas age-scale. The best age estimate (realization) is not necessarily the most frequent age estimate. Figure 4d shows an example from sample depth 621.28 m where there is a large difference between these two age estimates. Large differences can occur because the prior age estimate (i.e., the age estimate based only on prior ACPs) differs by a large amount from the “true” age of that sample and because the goodness-of-fit parameter considers the fit over the whole record. For depth 621.28 m most realizations resulted in an age estimate of 9200 years BP, similar to its prior age estimate of 9240 years BP, but the best realization estimated the age to be 9012 years BP. There are two possible reasons for this type of result: (1) this realization is the best because it managed to push the age of this sample by >200 years towards younger ages while not significantly changing the ages of nearby sections which already fit well or (2) no significant improvement in the goodness of fit was found by adjusting the age of this depth, and the goodness of fit was dictated by other sections of the record.

Uncertainty in gas age constraints is calculated as the root-mean-square error of FCP ages from the 1000 Monte Carlo realizations (Fig. 5d), although we acknowledge that this can overly simplify the empirical distribution of possible ages for a given depth (Fig. 4). Assessing uncertainty in this way integrates two types of error. The first is the ability to exactly match timing of two features. This uncertainty is determined by how well features are resolved in the records, measurement error, and the limited degrees of freedom in synchronization (i.e., the synchronization routine assumes constant accumulation during each subsections in a piecewise manner, although the true accumulation history is more likely to smoothly vary). The second type of error is analogous to deciding which feature in the reference record is the correct match. For example, the methane peak centered at 459.05 m (4675 years BP) is, in some realizations, matched to a peak at ∼4550 years BP in the WAIS Divide record, providing two distinct age possibilities (Fig. 4c).

Figure 5 RICE ice core age uncertainty. (a) Methane and (b) δ¹⁸O_atm records from the RICE (gray), WAIS Divide (Rhodes et al., 2015; Buizert et al., 2015b) (red), and NGRIP ice cores (Baumgartner et al., 2014) (green). (c) Sample spacing for methane (solid lines) and δ¹⁸O_atm (dashed lines) for the various cores. (d) Gas age uncertainty, relative to WD2014, for ages determined from a Monte Carlo analysis (solid triangles) and for extended gas age control points (open diamonds).

4.2 Extending the chronology with visually matched gas age control points: 30–83 ka

The more limited resolution of the RICE methane record below 719.3 m (30.66 ka) prevented use of the automated synchronization routine. Gas age control points were visually chosen between 719.30 m (30.66 ka) and 746.00 m (65.6 ka) and between 746.00 and 752.95 m (∼83 ka) (Table 1). Between 719.30 (30.66 ka) and 746.00 m (65.6 ka), ACPs were visually chosen by comparing the RICE discrete and CFA methane records and the δ¹⁸O_atm record to the WAIS Divide methane (Rhodes et al., 2015) and δ¹⁸O_atm records (Buizert et al., 2015b).

The WAIS Divide records and the WD2014 chronology end at 67.8 ka. For the deeper (older) section of the RICE ice core, between 746.00 and 752.95 m, we build target records from records of methane and δ¹⁸O_atm from the NGRIP ice core (Baumgartner et al., 2014; Landais et al., 2007, 2010). Using the NGRIP ice core is internally consistent with the WD2014 chronology, which is tied to a modified GICC05modelext age-scale between 30 and 67.8 ka (Buizert et al., 2015b). Because NGRIP is in the Northern Hemisphere the methane concentration is higher than that in the RICE record. To account for this interpolar difference, we scale the NGRIP target methane record (Baumgartner et al., 2014) to Antarctic values using an empirical least-squares fit between WAIS and NGRIP records between 55 and 67.8 ka.

Buizert et al. (2015b) found that the annual-layer-counted portion of the GICC05modelext chronology (0–60 ka) (Andersen et al., 2006; Svensson et al., 2008) is systematically younger than ages of corresponding features found in the U–Th absolute dated Hulu speleothem record (Reimer et al., 2013; Southon et al., 2012). A fit to Hulu ages was optimized by scaling the GICC05modelext ages linearly by 1.0063. This suggests that on average 6.3 out of every 1000 annual layers were not counted. For our NGRIP-based target records we adjust the NGRIP age-scale by adopting the scaling of Buizert et al. (2015b). Older ages in the GICC05modelext chronology are derived from the ss09sea Dansgaard–Johnsen model (Johnsen et al., 2001; NGRIP Community Members, 2004). We have added a constant 378 years (0.0063×60 000) for depths older than 60 ka in the target GICC05modelext ages to make this section continuous with the adjusted annual-layer-counted section.

Synchronization requires well-sampled identifiable features in both records. We estimate age uncertainty for visually matched ACPs (ACPs older than 29.9 ka) as the larger of the sample spacing of the two records being synchronized, following methods of Brook et al. (2005). Gas age uncertainty is plotted in Fig. 5d. The largest uncertainty is ∼1500–1700 years and occurs between 41.7 and 47.6 ka, where the age-scale is very compressed, and the RICE records are coarsely sampled.

4.3 Age control points in the disturbed ice: 746–753 m

Continuity of the RICE ice core appears to end at 746.00 m below surface, where a discontinuity of 20 ‰ is observed in the δD record (746.00–746.10 m) (Fig. 6b), and a 0.35 ‰ change is observed between sequential δ¹⁸O_atm samples (745.81 and 746.20 m) (Fig. 6c). The continuous chronology described in the previous section dates this ice to the end of DO-18 (∼64.6 ka). This age is supported by very negative δD values and a trend to more enriched values at shallower depths, which is consistent with warming trends observed in other Antarctic ice cores at this time period (Buizert et al., 2015b; EPICA Community Members, 2006; Petit et al., 1999; Parrenin et al., 2013) (Fig. 6f).

Figure 6 (a) Depth–age relationship from 740 to 756 m and evidence of an age reversal within the RICE ice core from (b) δD, (c) δ¹⁸O_atm, and (d) methane data. (b–d) Measurements are plotted against depth and color coded according to the age of the sample. Open circles represent samples for which an age could not be determined. (a–d) Red hatched bars represent discontinuities, where periods of climate history appear to be missing. Solid gray bars in panel (b) are measurement gaps in the δD record associated with large changes in δD. (e–h) The depth–age relationship and measurements are plotted against the age of the sample and color coded according to depth. WAIS Divide data (red), NGRIP (green), and the EDML δ¹⁸O_ice record (blue).

Figure 6 presents our best effort to date the ice in the disturbed section, between 746.00 and 752.95 m. The left panel of Fig. 6 shows the data as a function of depth, color-coded by age. The right panel shows the reconstructed time history of these variables, color-coded by depth. Below, we discuss how the ages are assigned to different depth sections (divided by vertical red lines in the left side of Fig. 6), starting with the shallowest portion. Dating folded ice remains challenging, and the solution we found may not be unique.

Immediately below the discontinuity, between 746.10 and 747.85 m, methane and δ¹⁸O_atm values best match our target records during DO-20 between 74.5 and 77.3 ka (Fig. 6), although the RICE methane record appears ∼30 ppb too high. Dating of this section to DO-20 implies either that 9.4 ka of climate history is absent from the RICE ice core because climate was not recorded or that 9.4 ka is compressed into ∼10 cm of ice. Extremely thin layers can be explained by ice flow or by an absence of accumulation, but the latter scenario would cause a collapse of the firn and δ¹⁵N-N₂ values approaching 0 ‰, which is not observed.

A cluster of samples between 747.85 and 750.46 m are ∼ 0.16 ‰ more enriched in δ¹⁸O_atm than was observed in the shallower section dated to DO-20 (Fig. 6c). In the NGRIP, EDML, and Siple Dome ice cores, a long-term trend in δ¹⁸O_atm towards more enriched values is observed from 80 to 65 ka (Capron et al., 2010; Severinghaus et al., 2009; Landais et al., 2007) (Fig. 6g). The enriched values between 748.34 and 750.46 m most likely indicate that this ice is younger than the adjacent shallower depths and that the stratigraphic order is reversed. This depth range is best matched to DO-19 between 71.8 and 74.3 ka (Fig. 6). Below 750.46 m, methane and δ¹⁸O_atm return to values that best match DO-20.

We note that two significant measurement gaps exist in the δD record, at 746.89–747.07 and at 750.00–750.25 m (vertical bars in Fig. 6b). A 12 ‰ shift in δD accompanies both of these gaps. The second measurement gap (750.46–750.56 m) is at nearly the same depth which separates the samples dating to DO-19 from the grouping immediately below dated to DO-20. The discontinuity of the δD record at these gaps may signify another hiatus in the climate record or highly contorted layers that are typically found around folds (Cunningham and Waddington, 1990; Alley et al., 1997; Thorsteinsson and Waddington, 2002; Waddington et al., 2001).

Below 750.46 meters, trends in methane and δ¹⁸O_atm indicate that the age of ice increases with depth until at least 753 m (∼83 ka). Age of this depth is constrained by a measured δ¹⁸O_atm value of −0.175 ‰. Such a depleted value is rare and only occurs during periods of high sea level and small ice sheets, the most recent time period prior to the Holocene being MIS-5a (80–85 ka) (Severinghaus et al., 2009; Capron et al., 2010; Petit et al., 1999; Landais et al., 2007; Kawamura et al., 2007). Negative values were observed at two other depths below 753 m (Fig. 2h), but the stratigraphic order of these depths is difficult to assess.

4.4 Firn modeling to determine Δage

Air transport in the firn causes an age difference between ice and air trapped in the ice, commonly referred to as Δage. Firn densification models are typically used to simulate past Δage (Schwander et al., 1997; Goujon et al., 2003). Input parameters (for example, temperature, accumulation rate, and surface density) are normally the largest sources of uncertainty. Using δ¹⁵N-N₂ as a proxy for past firn column thickness and assuming δD records past site temperature faithfully, firn densification models can be run in an inverse mode to estimate both past Δage and accumulation rates (Buizert et al., 2014, 2015b), which is the approach we employ here.

We use a dynamic version of the Herron–Langway model (Herron and Langway, 1980), which was also used for construction of the WD2014 chronology (Buizert et al., 2015b). The model simulates firn compaction rates and vertical heat diffusion and advection. The model domain is the full 763 m ice column at 0.5 m resolution; a time step of 1 year is used. The model simulates both gravitational enrichment of δ¹⁵N-N₂ and fractionation in the presence of thermal gradients. The model is forced with a temperature history derived from the δD record assuming a constant isotope sensitivity of 6 ‰ K⁻¹ (Brook et al., 2005). A smoothed (500 years LOESS filter) version of the δD data is shown in Fig. 3b and provided in the Supplement. The Δage model uses a heat diffusion model to determine firn temperatures. Differences between patterns of modeled temperature variation and that assumed from the smoothed δD record are minor (equivalent to Δage variations of decades, within model uncertainty). The assumed isotope-temperature sensitivity is similar to that observed for West Antarctica (Stenni et al., 2017; Cuffey et al., 2016; Masson-Delmotte et al., 2008) but somewhat higher than that obtained by comparing the RICE isotope record to ERA interim data (Bertler et al., 2018). The dependency of Δage on the assumed isotope sensitivity was explored in a model sensitivity test (Appendix C) and is incorporated into the Δage uncertainty estimates. The model further assumes a constant ice thickness, a constant 2 m convective zone similar to other high-accumulation dome sites with low mean wind speeds (Kawamura et al., 2006; Landais et al., 2006), surface firn density of 400 kg m⁻³ to match the modern surface firn density, and geothermal heat flux of 78 mW m⁻².

In a first iteration, we assume a prior ice age-scale for the temperature history by adding a constant 150 years to the gas chronology (Sect. 4.1 and 4.2). A new Δage solution is then calculated using the dynamic firn densification model. The Δage solution (and thereby the ice age-scale) is refined iteratively until it no longer changes appreciably (consistent within ∼1 year).

The climate at Roosevelt Island, with high accumulation and relatively warm temperatures, results in small Δage values and consequently relatively low absolute uncertainty in Δage compared to most Antarctic sites. Modern Δage (estimated at 60 m depth, the shallowest measurement of δ¹⁵N-N₂) is 146 years with a reconstructed lock-in depth (LID) of 48 m, consistent with the modern density profile (Bertler et al., 2018; Winstrup et al., 2019). Holocene Δage is small, ranging between 140 and 182 years. During the last glacial period, simulated Δage values fluctuate between ∼150 and 350 years.

Uncertainties in past Δage include the uncertainty in model inputs as well as the model itself. The uncertainty in the Herron–Langway model is conservatively estimated to be 20 %, based on differences between firn models (Lundin et al., 2017). We assessed the uncertainty due to model inputs using a steady-state Herron–Langway model that approximates the dynamic version but requires less computational time (Appendix C). In a sensitivity test, we randomly perturb the parameterizations of temperature and LID and assumptions of convective-zone thickness, surface firn density, and geothermal heat flux and recalculate the Δage solution (model parameters and their base values and ranges used in the sensitivity test are provided in Table C1). A total of 6000 iterations were run. The sensitivity tests include a wide range of temperature histories to account for the possibility that some variations in δD were caused by nonthermal effects such as variability in precipitation seasonality, moisture sources and pathways, and post-depositional vapor exchange. Model sensitivity is reported as the root-mean-square error of Δage calculations for each depth. Combined Δage uncertainty, provided in the included age-scale, is the root sum squares of the model uncertainty and the model sensitivity.

4.5 Comparison of gas-derived and layer counted chronologies: 0–2649 years BP

The gas-derived ice age-scale provides a chronology from 60 to 753 m depth (Fig. 2) that is independent of but overlaps the annual-layer-counted section (0–343.7 m; Winstrup et al., 2019). Figure 7 shows both chronologies and differences between them (positive values indicate that the annual layer counts are older than the gas-derived ages). The gas-derived ice age-scale agrees well with the annual-layer-counted age-scale, within 33 years, with similar trends in the implied annual layer thickness (Fig. 7c). Differences between chronologies can result from error in the synchronization of gas records, calculation of Δage for either the RICE or WAIS Divide ice cores, interpolation between ACPs, or annual layer counts in either core (because the WAIS Divide age-scale is used as a reference for the gas age-scale).

Figure 7 Comparison of the gas-derived (gray) and the annual-layer-counted ice age-scale (blue; Winstrup et al., 2019). Six absolute age markers (closed red circles) were identified, and an additional 67 volcanic horizons were cross correlated to volcanic horizons in the WAIS Divide core (plotted on WD2014 ice age-scale, open red circles; Winstrup et al., 2019). (b) Difference in ice age between the gas-derived ice age-scale and annual layer counts; positive values indicate that at the same depth the annual layer counts are older than in the gas-derived age-scale. Uncertainty estimates of the gas-derived ice age-scale (solid gray lines), the Δage estimate only (dashed gray lines), and the annual-layer-counted age-scale (blue dashed lines). (c) Interpretations of annual layer thickness from the gas-derived ice age-scale and annual layer interpretations.

The average age difference at depths of FCPs is −1 years (n=18, implied age from gas-derived ice age-scale being older than the layer counted chronology). Root mean square of the age difference is 17.3 years. Discrete points in the gas-derived ice age-scale can also be compared to the age of 67 volcanic peaks identified in the RICE ice core and correlated to peaks identified in the WAIS Divide ice core (open red circles in Fig. 7; Winstrup et al., 2019). Compared to these volcanic peaks, the root-mean-square ice age difference in the gas-derived ice age-scale is 13.6 years from their WD2014 ages. Good agreement between the two approaches gives confidence in the methodology used for the deeper section of the RICE17 chronology. The two largest differences occur at 89.72 m (243 years BP) and at 169.11 m (757 years BP) and are +30 and −33 years, respectively (Fig. 7b). These offsets are similar in magnitude to the individual uncertainties in calculating Δage or in synchronizing the gas records. The small age differences between the two ice chronologies also indicate that our approach to calculating uncertainty is likely conservative.

The RICE17 timescale transitions between the annual-layer-counted and gas-derived age-scales at 343.7 m, the deepest (oldest) FCP for which annual layers were identifiable (Sect. 4.1). The age difference at this depth is 3 years, with the gas-derived ice age implying an older age than the annual-layer-counted chronology. This age difference is much less than the respective age uncertainties of 45 and 111 years (2σ) for the annual layer counts and gas-derived ice age-scale, respectively. Because of the good agreement between the two ice age-scales, no correction is made for the 3-year offset, and the annual-layer-counted age is used.

5 Results and key observations from the RICE17 chronology and gas data sets

5.1 Implied accumulation history of Roosevelt Island

From 65 to 32 ka, the low variability in δ¹⁵N-N₂, ranging from 0.20 ‰ to 0.24 ‰, and the cooling trend observed for Antarctica for this time period suggest a decreasing trend in accumulation (Fig. 3).

While δ¹⁵N-N₂ values appear steady from 65 to 32 ka and during the Holocene (Fig. 3d), they exhibit large variability between 32 and 10 ka (Fig. 8a). δ¹⁵N-N₂ generally trends to heavier values beginning at 32 until 14.71 ka, indicating a growing firn column. An inflection point is observed at 25.3 ka, which is interpreted as an acceleration of firn thickening concluding in a peak of 0.293 ‰ at 21.8 ka. After a short depletion the steep trend resumed with a second δ¹⁵N-N₂ maximum of 0.326 ‰ reached at 15.7 ka. Following this maximum, δ¹⁵N-N₂ abruptly decreased from 0.308 ‰ to 0.220 ‰ at 14.71 ka, which corresponds with the beginning of the Antarctic Cold Reversal (ACR). At 12.38 ka, after the ACR, δ¹⁵N-N₂ partially recovers to pre-ACR values with an abrupt increase to 0.260 ‰.

Figure 8 Comparison of (a) δ¹⁵N-N₂, (b) CH₄, (c) annual layer thickness implied from depth–age scale, (d) lock-in depth (LID), and (e) accumulation reconstructions from the RICE ice core. Lock-in depth and accumulation are calculated with a dynamic Herron–Langway firn densification model.

Between 25.3 and 32 ka, accumulation is estimated to be ∼ 10 cm ice equivalent per year (Fig. 8e), and the increasing firn thickness is largely attributable to decreasing temperature. After 25.3 ka accumulation starts to increase, and by the first δ¹⁵N-N₂ maximum (21.8 ka), accumulation had increased to ∼17 cm ice equivalent per year. This feature is not apparent in other ice cores from the Ross Sea region, but those records tend to be difficult to interpret because of chronological uncertainties, such as in the case of Taylor Dome (Baggenstos et al., 2017), or because of unexplained jumps in δ¹⁵N-N₂, such as in the case of Siple Dome (Severinghaus et al., 2009). By the second maximum (15.7 ka), accumulation increased to ∼25 cm ice equivalent per year, similar to accumulation rates observed during the Holocene. An accumulation peak at the end of glacial terminations is consistent with evidence from trimlines in interior WAIS (Ackert et al., 2013) and is also observed in the accumulation histories from WAIS Divide and Siple Dome (WAIS Divide Project Members, 2013; Waddington et al., 2005) but to a smaller degree. The early accumulation peak at 21.8 ka is unique to the RICE ice core.

A large, rapid depletion in δ¹⁵N-N₂ is observed at 14.71 ka with low values lasting until 12.38 ka. Analysis of the lead–lag relationship between δ¹⁵N-N₂ and methane has been used to infer near-synchronous climate changes throughout the tropics and Northern Hemisphere (Rosen et al., 2014). These climate events are believed to originate in the Northern Hemisphere and propagate to the Southern Hemisphere (Blunier et al., 1998; Blunier and Brook, 2001; Buizert et al., 2015a, 2018). Curiously, the abrupt decrease in δ¹⁵N-N₂ at 14.71 ka (0.088 ‰, 683.70–681.80 m) is observed in samples deeper than the increase in methane marking the onset of the Bølling-Allerød (14.66–14.52 ka, 683.13–682.52 m), meaning that this climate event unambiguously precedes the Northern Hemisphere event (Fig. 9). Where as Rosen et al. (2014) specifically considered the thermal component of δ¹⁵N-N₂, the thermal component is considered to be negligible in Antarctic cores because air temperature only changes gradually. At Roosevelt Island, this period of low δ¹⁵N-N₂ values is interpreted to represent shallow firn thickness which in the firn model is caused by low snow accumulation (< 10 cm yr⁻¹; Fig. 8e). Low accumulation during this period is consistent with thin annual layers interpreted from the age–depth relationship (Fig. 3f), i.e., 0.3–0.6 cm yr⁻¹ annual layer thickness compared with 1.6–3.2 cm yr⁻¹ between 10.09 and 11.01 ka (642.75–661.07 m). Following the ACR, the modeled accumulation rate fully recovers to the ∼25 cm yr⁻¹ observed before the ACR.

Figure 9 Comparison of (a) methane and (b) δ¹⁵N-N₂ from the RICE ice core plotted on depth. Vertical blue bar highlights the onset of a large decrease in δ¹⁵N-N₂.

We propose three potential explanations for the thin annual layers and low accumulation rate observed during the ACR, with the last explanation being our preferred. (1) A large accumulation gradient is observed across the Roosevelt Island ice dome (Winstrup et al., 2019), which implies that a period of low accumulation at the RICE site may be the result of changes to the geometry of the Roosevelt Island ice dome. However, ice divide migration typically occurs over long timescales. (2) Accelerated ice flow may also cause thin annual layers and could potentially even affect layers within the firn. This flow could be the result of the ice streams which surrounded and buttressed Roosevelt Island ice dome suddenly being removed (Ackert et al., 1999, 2013; Halberstadt et al., 2016). The timing of the low-accumulation interval is similar to when dust records from the Taylor Dome ice core, a site located in the Transantarctic Mountains near the Ross Sea (Fig. 1), show an increase in dust diameter and a change in radiogenic isotope composition, which implies that the provenance of dust changed to a local source (Aarons et al., 2016). Aarons et al. (2016) interpreted this source as newly exposed land created by the withdrawal of Ross Ice Shelf. The timing of the low-accumulation period is also similar to some estimates of the timing of retreat of the WAIS based on sediment facies succession and radiocarbon dating (minimum date of 8.6 ka; McKay et al., 2016) and ice-sheet modeling (Golledge et al., 2014; McKay et al., 2016). Although annual layer thickness immediately above the ACR section is observed to be nearly 5 times as thick as the ACR section, this large change in thickness, if solely the result of changes in ice flow, would require an unrealistic change in the thickness of the dome or its ice flow which is not supported by TAC measurements or ice-flow modeling.

(3) Our preferred explanation is that an interval of low accumulation between 12.38 and 14.71 ka resulted in a shallow firn structure and depleted δ¹⁵N-N₂ values. In recent times, moisture arriving at Roosevelt Island is frequently related to enhanced cyclonic airflow in the Ross Sea and a strong Amundsen Sea Low (Tuohy et al., 2015; Emanuelsson et al., 2018). This period of low accumulation may indicate a changed atmospheric structure in the South Pacific where southward airflow is blocked by a persistent low-pressure zone north of the Ross Sea such as observed in more recent periods in the accumulation record from RICE (Bertler et al., 2018; Emanuelsson et al., 2015) and potentially related to the past ice shelf extent. Such a pronounced minimum in accumulation is not observed in the Siple Dome ice core, where annual layers are observed to thicken during the ACR (Brook et al., 2005; Waddington et al., 2005). If nonthermal effects influenced the RICE δD record, which is interpreted as temperature in our firn model, then the magnitude of accumulation change during this period may not be as large as reconstructed from the firn model. While not currently available, measurements of d_excess, dust particle size distributions, and dust geochemistry may be helpful in explaining the temperature and accumulation history of RICE.

5.2 Implications for climate and ice sheet history

Early reconstructions of the Ross Ice sheet during the LGM were based on glacial geological constraints from the western margin of the embayment. Denton and Hughes (2002) describe a maximum scenario in which thick ice in the Ross Embayment overrode both Siple Dome and Roosevelt Island. However, more recent observations and model experiments indicate a “fast and thin” scenario in which Siple Dome was not overrun by the interior ice sheet during the LGM (Waddington et al., 2005; Price et al., 2007), and although the Ross ice streams likely slowed down during the LGM, they remained active, maintaining a low elevation profile of the ice sheet in the Ross Sea (Parizek and Alley, 2004).

Our results further support the fast and thin scenario and add a key new constraint on ice thickness and thinning in the eastern Ross Sea. Specifically, our results suggest that ice deposited on Roosevelt Island originated as accumulation local to the drilling site, which may not be true if WAIS was thick during the LGM and overrode Roosevelt Island. If remnants of WAIS were stranded on Roosevelt Island, this would likely result in a hiatus in the gas and ice chronology, in values of δD or TAC characteristic of more continental precipitation or much higher elevations, or in discontinuities in the δ¹⁵N-N₂ record, indicating a much different firn structure.

While the RICE ice core chronology does exhibit an abrupt shift in δ¹⁵N-N₂ and TAC during the ACR (Fig. 3), no discontinuity in δD was observed, meaning that it is unlikely that any of this ice originated as part of WAIS. At least one age discontinuity (at 64.6 ka) as well as an age reversal was observed deeper in the core, but these depths do not coincide with the timing of the retreat of WAIS in the Ross Sea (McKay et al., 2016). The largest discontinuity, at 746.00–746.10 m, is accompanied a 20 ‰ change in δD. Methane and δ¹⁸O_atm indicate that this discontinuity represents a 9.4 ka age gap. Over the same age range, the EDML ice core records a similar change in δD (Fig. 6f), indicating that this change in δD is explained by Antarctic climate patterns alone and without invoking large changes in ice sheet configuration. A possible second discontinuity is observed at 747.00 m depth at a small gap in measurements of 2.7 cm. Ice immediately below 747.00 m is interpreted as being warmer, meaning that this is probably not derived from somewhere upstream in WAIS or from a higher elevation (Fig. 6b). Additionally, no dramatic or sudden change in δ¹⁵N-N₂ or TAC was observed in association with either of these discontinuities (Fig. 6). We conclude that it is highly unlikely that the accumulation site of the RICE ice core changed during the deglaciation and that Roosevelt Island ice dome probably remained independent of an advanced WAIS during the LGM. A similar conclusion was reached about Siple Dome during the LGM (Waddington et al., 2005; Nereson et al., 1998; Price et al., 2007; Parizek and Alley, 2004).

Geomorphological features on the Ross Sea bed do provide evidence of grounded ice north of Roosevelt Island during the LGM (Shipp et al., 1999; Anderson et al., 1984, 1992, 2014; Halberstadt et al., 2016). If these features were formed by an extended WAIS, it would imply that ice flowed around Roosevelt Island and Siple Dome and therefore was limited in thickness. These conditions would indicate that ice streams were active throughout the last glacial period. Alternatively, these geomorphic features may be the result of ice from a different origin. Price et al. (2007) proposed that during the LGM, an ice dome may have existed on Mary Byrd Land. In this scenario, thick, grounded ice could exist north of Roosevelt Island without flowing over or around the Roosevelt Island sea rise. The RICE records can not distinguish between these scenarios.

6 Conclusions

We present the RICE17 gas and ice chronologies for the RICE ice core. These timescales date the gas and ice records for the last ∼83 ka. Between 0 and 30 ka an automated synchronization routine is used to identify gas age control points that best match the RICE methane and δ¹⁸O_atm records to the respective records from the WAIS Divide ice core (WAIS Divide Project Members, 2013; Buizert et al., 2015b). This technique requires few prior constraints, accommodates simultaneous synchronization of multiple parameters, and allows assessment of age uncertainty. Unique in this approach is the use of centennial-scale variability in methane for chronostratigraphic matching for ages older than the last ∼2000 years. Synchronization between ice cores for the time period between 30 and 83 ka (719–753 m) was accomplished by manually choosing tie points. Below 753 m the ice could not be dated with the currently available data. The RICE17 ice age-scale is based on annual layer counts between −62 and 2649 years BP (0–343.7 m), and for depths below 343.7 m, a separate ice age-scale was derived from the synchronized gas age-scale by adding Δage estimated from a firn model.

A key contribution from the development of the RICE17 age-scale is evidence of active ice streams in the eastern Ross Sea during the last glacial cycle. This is supported by the continuous age-scale and continuous records of climate from RICE which imply that the Roosevelt Island ice dome remained stable and independent of WAIS for at least the last 64.6 ka and likely for the last 83 ka.

The RICE ice core provides records of climate, with precise dating, for a scarcely sampled region of Antarctica. Future work will investigate regional climate of the eastern Ross Sea in comparison to climate records from other sites to better understand spatial patterns around Antarctica, such as during the ACR when Roosevelt Island experienced an interval of particularly low accumulation, and to study the glacial retreat of WAIS in the Ross Sea.

Appendix A: Methane measurements

The methodology used at Oregon State University for measuring methane concentration in ice cores was described in Grachev et al. (2009) and Mitchell et al. (2013). Briefly, 40–60 g ice samples are trimmed and placed in glass flasks. The glass flasks are then attached to a high vacuum line and immersed in a chilled ethanol bath set to −63 ^∘C to keep the samples frozen. Since Mitchell et al. (2013), insulation has been added around the ethanol bath and above where the flasks are mounted. The added insulation reduced the temperature and water vapor content of gas in the headspace of the flasks and decreased sample-to-sample temperature variability. Both can affect methane measurements by changing the pressure reading or the retention time of methane in the gas chromatograph (GC) column. We also made efforts to more carefully regulate the amount of ethanol in the chilled bath and temperature of the hot water bath during melting. These steps improved stability of gas extraction and of measurements between sampling days.

After laboratory air has been evacuated, the flask valves are closed, and the samples are melted in a hot water bath for 30 min to liberate air. Samples are refrozen by immersing the flasks in the cold ethanol bath. Once the samples are completely refrozen and the sample flask temperature has stabilized (approximately 1 h after refreezing begins), methane measurement is performed by expanding sample air from the flask headspace into a gas chromatograph (GC).

Calibration measurements are made on an internal standard which is referenced to several compressed air standards externally calibrated on the NOAA scale (Dlugokencky et al., 2005). Calibration runs of systematically varying pressures are made both before and after samples are measured. The calibration curve is a least-squares linear regression between pressure and peak area for both sets of calibration runs described by a slope m_cal and intercept b_cal.

Methane concentration (C) was previously calculated by comparing the sample pressure (P_meas) and the peak area (PA_meas) from GC analysis to the predicted peak area of the standard gas of equal pressure.

$$\begin{array}{}\text{(A1)}& C_{\mathrm{raw}}=C_{\mathrm{standard}}\cdot {\displaystyle \frac{{\mathrm{PA}}_{\mathrm{meas}}-b_{\mathrm{cal}}}{P_{\mathrm{meas}}\cdot m_{\mathrm{cal}}}}\end{array}$$

Four measurements can be made on each ice sample. If the two sets of calibration runs are taken individually, interpretation of sample concentrations typically differ by <3 ppb.

In this calculation, the sample pressure is used to quantify the number of moles of air in the sample assuming that the sample temperature remains constant and equal to standard air temperature during calibration runs. However, in a series of dry blank experiments it was determined that sample gas temperature is cooler than the standard gas temperature. We estimate the relative temperature difference to be 0.42 % of the standard air temperature and thus the sample pressure needs to be corrected upwards by that amount. We call this the GC-thermal effect. We now apply the correction to the previous concentration calculation.

$$\begin{array}{}\text{(A2)}& C_{\mathrm{raw}}=C_{\mathrm{standard}}\cdot {\displaystyle \frac{{\mathrm{PA}}_{\mathrm{meas}}-b_{\mathrm{cal}}}{P_{\mathrm{meas}}\cdot \mathrm{1.0042}\cdot m_{\mathrm{cal}}}}\end{array}$$

Results from typical wet blank samples, in which we add standard air over bubble-free ice made from Milli-Q ultrapure water, had shown that these samples are typically 2–3 ppb enriched compared to the standard concentration using the old method for calculating concentration. This agrees with the predicted magnitude of enrichment from the GC-thermal effect of 2.1 ppb.

Several wet blank samples are measured each day and are interspersed between ice core samples. The offset between the measured concentration of the wet blank and the known concentration of the standard gas added is called the “blank” correction. This represents any effects during the sample analysis process which may alter the measured concentration. We bin the wet blank results over the time period for which samples were measured to establish a single blank correction value with an estimated blank correction uncertainty. RICE ice core samples were measured during two separate periods and have separate blank correction values.

Since the OSU analytical methods use a wet extraction technique to liberate ice core air, effects of gas solubility must be accounted for. Methane is more soluble than the major components of air and is preferentially dissolved during the extraction step. This leads to a decrease in the measured concentration compared to the true ice core concentration. Mitchell et al. (2013) empirically derived a solubility correction (S) which we repeated several times for the RICE samples. S is defined as the total amount of methane (gas + dissolved) divided by the amount in gas phase. A solubility value of 1.0165 was used for ice core samples and 1.0079 for bubble-free ice. We believe the difference in S for the different ice samples results from the differences between how blank ice and ice containing air behave during melting. Specifically,

• Bubble-free ice melts slowly in comparison to glacial ice, which sometimes melts rapidly and cracks violently. This, along with bubbles rising and breaching the meltwater surface, causes disturbances in the water–air interface and promotes exchange of CH₄ into the meltwater. This should lead to greater mixing and homogenization of air and water.

• Bubbles released into the meltwater will be at higher pressures than the overlaying air because of surface tension. The higher partial pressure of CH₄ in those bubbles, in comparison to the standard gas added over the bubble-free ice, will cause air to go into solution faster.

• Because glacial ice tends to be melted sooner than blank ice, a longer time period for liquid–gas exchange is available.

The blank and “solubility” corrections are applied in the following way:

$$\begin{array}{}\text{(A3)}& C_{\mathrm{corrected}}=C_{\mathrm{raw}}\cdot S_{\mathrm{sample}}-\left(C_{\mathrm{blank}}\cdot S_{\mathrm{blank}}-C_{\mathrm{standard}}\right).\end{array}$$

This formula differs from that used by Mitchell et al. (2013). In Mitchell et al. (2013) no solubility correction was applied to the bubble-free ice samples. The difference results in a constant offset of −7.4 ppb compared to Mitchell et al. (2013).

Appendix B: Calibration of RICE CFA methane

The RICE CFA methane record is qualitatively used for synchronization to the WAIS Divide methane record, and thus careful calibration was not required. Regardless, we present an ad hoc calibration of the RICE CFA methane record based on comparison to the RICE discrete methane record measured at Oregon State University (OSU). The RICE CFA methane record was measured in multiple campaigns, and major adjustments to the analytical system occurred during both of those periods. Calibration of the RICE CFA record is done in a piecewise manner to reflect these changes.

Our calibration scheme accounts for instrument calibration, a concentration-based correction due to instrument sensitivity, and measurement drift, a sampling-time-dependent correction. For comparison between data sets, we subsample the CFA data at depths where discrete measurements were made. We first apply the concentration calibration by regressing the subsampled CFA methane concentration against the discrete measurement. Drift was accounted for by a second regression comparing the residual of the concentration calibration against either the measurement time or sample depth, which ever provided the better correlation to the discrete data set. Drift was only a significant factor between 500 and 680 m depth. Calibration values are given in Table B1.

Table B1 Calibration of the RICE CFA methane data set is performed by comparison to the RICE discrete methane data set with corrections for instrument sensitivity and instrument drift. Calibration is done for different segments of the core. Quality of the fit is described by the R² statistic comparing calibrated values of the subsampled CFA methane record to discrete measurements.

Uncertainty in the relationship may be caused by measurement uncertainty, which is relatively small, or the uncertainty in the depth registration of the continuous measurements. Uncertainty in the depth registration is relatively unimportant in the top ∼670 m of core. In this section annual layers are relatively thick, and methane variability is relatively low, which results in only minor uncertainty in methane concentration from an error in depth. However, for the deeper section of core, methane concentrations vary rapidly with depth, and small errors in the depth registry represent large differences in methane concentration. Because of errors in the depth registry and the sensitivity of the inferred methane concentration, we restrict our calibration scheme to the last ∼39.66 ka (725.63 m) ending after the inclusion of GIS-8.

Appendix C: Steady-state Herron–Langway sensitivity test

The gas age–ice age offset (Δage) for the RICE ice core was estimated with a dynamic Herron–Langway model (Buizert et al., 2015b). The model is constrained by measurements of δ¹⁵N-N₂ (a proxy for firn thickness) and δD (a proxy for past temperature). The dynamic model also assumes a convective zone of 2 m, surface firn density of 400 kg m⁻³, and geothermal heat flux of 78 mW m⁻².

A steady-state version of the Herron–Langway model mimics the dynamic version following a similar iterative methodology. A constant Δage estimate is assumed as a prior to assume an ice age-scale for the temperature history. The model then calculates a new Δage solution from the temperature and firn thickness histories. This process is repeated until iterations no longer change appreciably.

While the dynamic version is used to establish the Δage history for RICE, the steady-state version has the advantage of being computationally faster and is used in a sensitivity test. In this test, we vary prior assumptions about the isotopic temperature sensitivity used to infer past temperature, the convective zone thickness, surface firn density, and assumptions about geothermal induced temperature gradient in the firn. In a Monte Carlo approach, each parameter is given a range of acceptable values from which the steady-state Herron–Langway model calculates the Δage history (Table C1). This is repeated for 6000 realizations. Realizations are rejected if any of the following points are true:

Table C1 Model Parameters used for steady-state Herron–Langway model.

• The modeled modern ice age at lock-in depth (LID) is more than 25 years different than the annual-layer-counted age of that depth (48.57 m, 89 years BP).

• The modeled modern accumulation is less than 0.15 or greater than 0.35 m ice per year.

• The isotopic sensitivity is less than 3.2 or greater than 9.6 ‰ K⁻¹, similar to the range of values observed for West Antarctica (Masson-Delmotte et al., 2008; Cuffey et al., 2016; Stenni et al., 2017). This parameter is randomly chosen from a normal distribution at the beginning of each realization and can fall outside of this range.

• The minimum temperature is less than −60^∘C.

• The maximum estimated Δage is more than 1000 years.

• The maximum modeled accumulation exceeds 1.0 m ice per year or minimum accumulation is negative.

Δage is calculated for each depth that δ¹⁵N-N₂ was measured. Uncertainty in Δage is assumed to be the root-mean-square error of accepted realizations and is estimated for each sample depth.