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

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 2Data 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.

Download

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 3Data 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.

Download

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.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.

Download

Table 1Prior 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.

Download Print Version | Download XLSX

“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 2Model Parameters used for optimized correlation routine. Code for routine can be found in the Supplement.

Download Print Version | Download XLSX

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 5RICE 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).

Download

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).

Download

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 7Comparison 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.

Download

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.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 8Comparison 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.

Download

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 9Comparison 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₂.

Download

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.