2.1. South Pole site specifics

The main coring of the South Pole Ice Coring Project (SPICEcore) took place over two Antarctic summer seasons between 2014 and 2016. The drilling site, located at 89.99° S, 98.16° W and ∼2 km from the primary South Pole Station, is unique in that the mean annual temperature at the drilling location (−49°C) is more typical of East Antarctic sites, while the mean accumulation rate (∼7.5 cm w.e. a⁻¹) is more typical of West Antarctic sites (Casey and others, Reference Casey, Fudge, Neumann, Steig, Cavitte and Blankenship2014). The site has an estimated ice thickness of 2700 m, an ice velocity of 10 m a⁻¹ along 40°W, is situated at ∼2835 m above sea level, and has an estimated depth to pore close-off of ∼120 m (Battle and others, Reference Battle1996, Casey and others, Reference Casey, Fudge, Neumann, Steig, Cavitte and Blankenship2014; Sneed and others, Reference Sneed, Mayewski and Dixon2011). The primary recovered core, SPC14, was drilled to a depth of 1751 m, with the ice dated to ∼54 300 years bp (Winski and others, Reference Winski2019). Calculated gas-age/ice-age differences for the core were as high as 2500 years, with gas-ages at the bottom of the core of ∼52 500 years bp (Epifanio and others, Reference Epifanio2020).

Several studies have been published recently reconstructing glaciological, chemical and temperature histories for the South Pole site using the primary SPC14 core. These include analyses of accumulation (Fudge and others, Reference Fudge2016; Kahle and others, Reference Kahle2021), advection and ice-provenance trends (Lilien and others, Reference Lilien2018; Fudge and others, Reference Fudge2020), oxygen isotope trends (Markle and others, Reference Markle, Steig, Roe, Winckler and McConnell2018; Steig and others, Reference Steig2021; Markle and Steig, Reference Markle and Steig2022) and trapped gas trends (Epifanio and others, Reference Epifanio2020; Morgan and others, Reference Morgan, Buizert, Fudge, Kawamura, Severinghaus and Trudinger2022; Epifanio and others, Reference Epifanio2023). The temperature reconstructions published thus far for the South Pole site have used several proxies, including the well-known water stable isotope thermometer (Markle and Steig, Reference Markle and Steig2022), isotope diffusion lengths (Kahle and others, Reference Kahle2021) and the history of firn densification at the site (Buizert and others, Reference Buizert2021a). These reconstructions for South Pole have all revealed a general temperature increase across the recent glacial transition of between ∼7°C and ∼10°C.

An investigation of the visual stratigraphy was also carried out for the site (Fegyveresi and others, Reference Fegyveresi, Fudge, Winski, Ferris and Alley2019a). Seasonal coarse-grained and depth-hoar layers were used to identify the annual layers (also see Gow, Reference Gow1965), and the results reveal accumulation trends that could indicate past climatological or glaciological changes at the site. Specifically, one 1600 year interval from 6700 to 5100 years bp (years bf 1950) experienced a higher-than-average accumulation rate (∼8.1 cm w.e. a⁻¹), while one 700 year interval from 3100 to 2400 years bp (years bf 1950) experienced a lower-than-average accumulation rate (∼6.4 cm w.e. a⁻¹); several analyses have suggested these temporal patterns in the core are related to upstream spatial patterns in accumulation and redistribution of snow (Winski and others, Reference Winski2021, Morgan and others, Reference Morgan, Buizert, Fudge, Kawamura, Severinghaus and Trudinger2022).

Over 1900 individual wind or iced crusts were also observed and documented, as well as one visible tephra layer at a depth of 306.6 m in the core (∼3560 years bp and likely tied to an eruption of Candlemas Island in the South Sandwich Islands; Palais and others, Reference Palais, Kyle, Mosley‐Thompson and Thomas1987). These data contributed to the accumulation-rate reconstruction that is needed in the BND temperature-reconstruction technique, and demonstrated, not surprisingly, that the core is ordinary bubbly ice suitable for that reconstruction.

2.2. BND paleothermometry

Even with the demonstrated success of known ice-core paleoclimate proxies (such as δ¹⁸O_ice), it is essential that any such proxy or its associated reconstruction technique or model be continually tested, refined, and improved to ensure that they are spatially and temporally robust across sites with differing paleo-histories and conditions. It is also valuable to develop complementary paleoclimate reconstructions using methods reliant on entirely independent variables, so as to better corroborate those reconstructions and affirm that they represent the most accurate temperature histories. A recent study showed that borehole thermometry in interior East Antarctica yields Last Glacial Maximum (LGM) temperature estimates significantly warmer than those estimated from the traditional use of water stable isotope ratios (Buizert and others, Reference Buizert2021a), providing a need for independently calibrated temperature proxies to resolve this discrepancy.

Based on the modeling in Alley and Fitzpatrick (Reference Alley and Fitzpatrick1999), Spencer and others (Reference Spencer, Alley and Fitzpatrick2006) developed a technique that allows for the reconstruction of either paleotemperatures or paleo-accumulation rates by fitting a semiempirical steady-state model to measured number-densities of bubbles in ice cores. Polar firn densification is primarily driven by temperature and accumulation rate (Gow, Reference Gow1968); however, during burial following pore close-off, ice grain sizes quickly become unreliable indicators of firnification conditions due to various processes that depend on other factors such as ice deformation rate (Spencer and others, Reference Spencer, Alley and Fitzpatrick2006). Given that the geometry of firn at pore-close off is essentially self-similar (Gow, Reference Gow1968), the integrated effects of temperature and accumulation rate on density and grain growth are therefore recorded in the number-density of bubbles formed as the firn is transformed into ice at the pore close-off depth (Spencer and others, Reference Spencer, Alley and Creyts2001, Reference Spencer, Alley and Fitzpatrick2006). Furthermore, due to slow gas diffusion between ice-core bubbles, BNDs in ice can reliably preserve information about firnification conditions over longer periods (Alley and Fitzpatrick, Reference Alley and Fitzpatrick1999).

Conserved number-density then allows for the use of the firn-densification and grain-growth models to invert for either the firnification paleotemperature or paleo-accumulation rate, provided the other parameter is known. This method was updated and applied successfully using ice-core data from the West Antarctic Ice Sheet (WAIS) Divide site, revealing an overall ∼1.7°C Holocene cooling—with notable smaller climatic trends contained therein—that also agreed with δ¹⁸O reconstructions (Fegyveresi and others, Reference Fegyveresi2011, Reference Fegyveresi2016).

As noted, the South Pole site experiences a unique combination of environmental variables for an Antarctic site, most notably its relatively high mean-annual accumulation rates for its relatively low mean-annual temperature (Casey and others, Reference Casey, Fudge, Neumann, Steig, Cavitte and Blankenship2014). Both of these site conditions favor slower grain growth with depth, and therefore smaller mean grain and bubble sizes, as well as higher BND values (Gow, Reference Gow1968). Combined with a higher-resolution annual stratigraphy as compared to other East Antarctic sites, it places South Pole into a ‘sweet spot’ with respect to BND, allowing for both a greater depth of resolvable annual layers, as well as the potential reconstruction of temperature histories back through the most-recent glacial termination into the Pleistocene (the Spencer model has not yet been successfully applied to any ice core containing bubbly ice older than Holocene-aged). It is worth noting that the South Pole site is located in a flank-flow setting, and therefore the climatic record from deeper in the core represents conditions from farther upglacier, introducing the need for ice-flow corrections to answer certain paleoclimatic questions, with the associated uncertainties in ice flow (Fudge and others, Reference Fudge2016; Lilien and others, Reference Lilien2018) .

During field recovery of the SPC14 core, it was also observed that the core exhibited very minimal brittle-ice behavior, even at brittle-ice zone (BIZ) depths traditionally associated with extreme instability, fracturing and spalling (Souney and others, Reference Souney2021; Fegyveresi and others, Reference Fegyveresi, Fudge, Winski, Ferris and Alley2019a). Later, inspection of thin sections prepared from the core also revealed a nearly complete lack of any chessboard subgrain structures, likely arising in large part from the small grain sizes which in turn are controlled by the high accumulation rate and low temperature of the site. Recent work by Fitzpatrick and others (Reference Fitzpatrick, Wilen, Voigt, Alley and Fegyveresi2024) and Barnett and others (Reference Barnett, Fegyveresi, Alley, Smith, Fitzpatrick and Voigt2023) suggests that brittle behavior involves fracturing along the chessboard subgrain boundaries. This diminished brittle-ice behavior has implications for site-selection for future ice cores. It also was extremely favorable for more continuous and accurate BND measurements, as bubble-section preparation was less affected by cracking during the ice cutting and microtoming processes, resulting in more-pristine (and defect-free) imagery.

Because of the higher BND counts, lack of overall brittleness and smaller grains, the South Pole core is excellent for both applying and refining the Spencer model, and for testing new measurement techniques and methods. As noted by Fegyveresi and others, (Reference Fegyveresi2011, Reference Fegyveresi2016), measurement of BND on bubble-sections introduces small uncertainties due to the two-dimensional (2-D) stereological projection of 3-D bubbles, which might be avoided with 3-D measurement techniques that capture true volumetric data and imagery, such as X-ray microtomography.

Over the past two decades, micro-computed tomography, or simply micro-CT, has been shown to be an effective tool for measuring 3-D properties of glacier ice without significant degradation to the sample (e.g. Hagenmuller and others, Reference Hagenmuller, Chambon, Lesaffre, Flin and Naaim2013; Dadic and others, Reference Dadic2019; Bagherzadeh and others, Reference Bagherzadeh, Freitag, Frese and Wilhelms2023). Micro-CT technology generates 3-D images composed of an integrated collection of 2-D projections, or ‘slices’, of the target specimen. At high enough resolution, this technique can therefore generate not just true 3-D counts of trapped bubbles in a given volume, but also volumetric size estimates of those bubbles. Here, we prepared 11 unique micro-CT samples at 100 m intervals (and from equivalent depths to those of traditional bubble-section sampling) as a way to test the efficacy and accuracy of the technology for measuring BND and mean bubble sizes. Our primary objective was to determine if this technique could be used as a viable alternative for ascertaining bubble counts, and ultimately, BND. We show that this new technique is viable and indeed nondestructive of the prepared samples, and can greatly reduce sample preparation and imaging processing time for future ice-core measurements.

3.1. Bubble number-density

To obtain both the ice-core bubble-sections and micro-CT sections, samples were first cut from the appropriate core depth sections in the National Science Foundation Ice Core Facility’s −26°C ice preparation room. Following Fegyveresi and others (Reference Fegyveresi2011, Reference Fegyveresi2016), we used a band saw to first cut thick bubble sections ∼10 cm long by ∼6 cm wide by 5 mm thick, at 20 m depth intervals ranging from 160 to 1200 m depth in the core. One face of each sample was smoothed with a sledge-type microtome and then affixed to a glass slide using a combination of silicone oil and tinkered water-ice. These samples were then microtomed to ∼1.5 mm overall thickness and recorded with several overlapping high-resolution images using the on-site digital camera stage. Final imagery was stitched using Adobe Photoshop software to produce the high-resolution bubble-sections (e.g. Fig. 1). Estimates of BND were calculated—with appropriate bubble-cut corrections—following Fegyveresi and others (Reference Fegyveresi2011, Reference Fegyveresi2016), and using the QIA-64 (formerly Fovea Pro) image analysis plug-in for Adobe Photoshop.

Figure 1. An unedited bubble-section image from the 160 m vertically oriented sample (just below pore close-off). Image shown with scale and stratigraphic core orientation.

Micro-CT samples were cut at 100 m intervals ranging from 200 to 1200 m depth. These samples measured 2.5 cm square, by 8 cm long and were split in half for replicate measurements. All micro-CT measurements were made on site in a − 10°C cold room at the Cold Regions Research and Engineering Lab in Hanover, NH. We used a Bruker SkyScan 1173 desktop micro-CT scanner with a Hamamatsu 130/300 tungsten X-ray source and flat-panel sensor camera detector, with samples mounted in standard plastic measurement containers (see, e.g. Fig. 2). The X-ray source of the SkyScan 1173 produces a fixed conical, polychromatic beam that operates a distortion-free X-ray detector with a source spot size of <5 µm, a spatial resolution detectability of <4 µm, and contrast resolution of <7 µm. For our measurements, we set the maximum accelerating voltage of the X-ray beam at 40 kV, with a current of 200 μA. Samples were rotated 180° in 0.35° steps, with five-frame averaged attenuation images captured at each step using a camera exposure of 400 ms. We used a 2 × 2 binning protocol to create final X-ray radiographs of 1120 × 1120 pixels. Reconstructions and analyses of these radiographs were carried out use the native Bruker software package NRECON, which uses a modified Feldkamp cone-beam algorithm to produce a vertical stack of gray-scale cross-section images. As part of the image post-processing, we performed ring artifact reduction, post-alignment and beam hardening corrections, and a two-pixel Gaussian kernel smoothing to reduce noise. The resulting images had a spatial resolution of 15 μm per voxel (e.g. Fig. 3). We selected a smaller internal volume of interest for analyses, measuring ∼10 000 mm³, in order to eliminate sample edge effects and gaps. We performed 3-D analyses on the resulting segmented images using Skycan CTAn software. Standard binary thresholding and 5-voxel despeckling algorithms were first applied to reduce image noise and eliminate non-bubble defects. Final returned data included all measured feature (bubble) counts and sizes. Overall measurement errors reflect a combination of analytical instrument uncertainty and the one sigma (1σ) standard of replicate samples.

Figure 2. Bruker SkyScan 1173 X-ray micro-CT instrument located in a −10°C cold room at the Cold Regions Research and Engineering Laboratory. Scanner shown with a sample mounted in plastic specimen container on the measurement stage.

Figure 3. A volumetric 3-D rendering of micro-CT imagery for the sample taken from 300 m depth. Each small feature represents a bubble measured within the sample volume of interest.

BND is the number of bubble centers in a given volume of ice, but because of the finite size of bubbles, parts of some bubbles with centers outside of a sample will be observed in a sample. Standard correction methods have been developed for these ‘cut bubbles’, and were followed here (see Saltykov, Reference Saltykov1976; Martinerie and others, Reference Martinerie, Lipenkov and Raynaud1990; Fegyveresi and others, Reference Fegyveresi2011). We also followed Lipenkov (Reference Lipenkov and Hondoh2000) and Ueltzhoffer and others (Reference Ueltzhöffer2010) in identifying and eliminating any ‘microbubbles’ from our bubble counts. In ice from especially cold sites such as South Pole, microbubbles can account for up to 20% of total observed bubbles (Lipenkov, Reference Lipenkov and Hondoh2000), although in our measured samples, they were never >5%. Consequently, the highest concentration of microbubbles we eliminated in any of our samples was 5% of the total bubble count, and errors in identification are a small fraction of that. Similar to WAIS Divide (Fegyveresi and others, Reference Fegyveresi2016), for the purpose of using BND as a climate indicator, this small concentration was negligible.

We measured BND using the bubble-sections and micro-CT samples, with typical values ranging between 800 and 900 bubbles cm⁻³ over the clathrate-free portion of our sampling interval—with a maximum value of 938 bubbles cm⁻³ measured at depth of 500 m (Fig. 4). These number-densities are notably higher than published values from other ice-core sites (see, e.g. Spencer and others, Reference Spencer, Alley and Fitzpatrick2006; Fegyveresi and others, Reference Fegyveresi2011), reflecting the combined effects of the relative colder mean temperatures, and higher mean accumulation rates at South Pole. The onset of clathrates did begin to notably affect our data below ∼1100 m, so for the purposes of our paleoclimate reconstruction, modeled temperatures from below this depth were considered unreliable. The glacial–interglacial transition, occurring between ∼1100 and 800 m in depth, is not immediately obvious in the raw BND data. This is because the colder temperatures and lower accumulation rates during glacial periods have opposite effects on BND.

Figure 4. Measured bubble number-densities plotted against depth. Error bars represent estimated mean bubble count error across all samples. The clathrate-ice zone is shown in gray.

3.2. Accumulation-rate history

As noted above, reconstruction of paleotemperatures from BND data requires paleo-accumulation rates. Here, we describe our accumulation history after a short review of earlier work illustrating important characteristics of the site.

Beginning in 1965, several shallow firn and ice cores were drilled near South Pole Station (e.g. Gow, Reference Gow1965; Giovinetto and Schwerdtfeger, Reference Giovinetto and Schwerdtfeger1965; Kuivinen, Reference Kuivinen1983; Mosley-Thompson and others, Reference Mosley-Thompson, Paskievitch, Gow and Thompson1999), each resulting in slightly different modern mean accumulation-rate estimates, ranging between 5.0 and 8.5 cm w.e. a⁻¹ (see also van der Veen and others, Reference Van der Veen, Mosley-Thompson, Gow and Mark1999; Casey and others, Reference Casey, Fudge, Neumann, Steig, Cavitte and Blankenship2014). During the 2002–03 field season, the International Trans-Antarctic Scientific Expedition (ITASE) drilled a 207 m core ∼8 km from the South Pole Station as part of their array of Antarctic cores. Visual stratigraphy was measured on this core in the field, yielding an accumulation rate of ∼7.6 cm w.e. a⁻¹ when averaged over the 207 m (A. Gow, pers. comm., 2015; Sneed and others, Reference Sneed, Mayewski and Dixon2011). The most recent estimate of modern mean-accumulation rate derived from the SPC14 core is ∼7.4 cm w.e. a⁻¹, ranging between ∼7.0 and ∼9.0 cm w.e. a⁻¹ through the Holocene (Fegyveresi and others, Reference Fegyveresi, Fudge, Winski, Ferris and Alley2019a; Winski and others, Reference Winski2019).

Because the South Pole site is in a region of flank flow (not on an ice divide), with notable uncertainties in flow velocities over the history of accumulation of the ice in the core, the origin of measured accumulation can be somewhat complicated (Casey and others, Reference Casey, Fudge, Neumann, Steig, Cavitte and Blankenship2014; Lilien and others, Reference Lilien2018). Estimated regional ice-flow velocities for South Pole range between 9.6 and ∼10.1 m a⁻¹ (Hamilton Reference Hamilton2004; Bamber and others Reference Bamber, Vaughan and Joughin2000), and therefore the provenance of any measured ice from deep within recovered cores at the site, is likely hundreds of kilometers upstream and sourced from a very large total catchment area.

Deeper layers in the SPC14 ice core were deposited upstream of the current coring location; therefore, the effective accumulation rate (A) in the core reflects both past climatic changes and upstream patterns of snow accumulation (Lilien and others, Reference Lilien2018; Fudge and others, Reference Fudge2020). Variations in the latter appear well-correlated with the surface curvature, or second derivative of the surface elevation, along the flow path, suggesting a role of surface snow redistribution by katabatic winds (Morgan and others, Reference Morgan, Buizert, Fudge, Kawamura, Severinghaus and Trudinger2022). The SPC14 (A) history has been reconstructed in various ways. Winski and others (Reference Winski2019) reconstructed Holocene accumulation rates based on annual layer thickness and a simple 1D Nye ice flow model (Nye, Reference Nye1963). The advantage of this method is that it allows high temporal resolution; however, it relies on a simplified and somewhat poorly constrained ice-flow model. Kahle and others (Reference Kahle2021) used a statistical inverse approach to simultaneously reconstruct a number of glacio-climatological parameters (A, T, thinning function) using empirical reconstructions of the water isotope diffusion length and the gas age-ice age difference (Δage). The advantage of this approach is that it covers the full length of the ice core; however, the method is not optimized for reconstructing past A, and it does not incorporate firn thickness information available from the isotopic composition of N₂ (δ¹⁵N-N₂). Here we use an empirical method to reconstruct A based on its fundamental relationship to firn thickness and firn age (Buizert, Reference Buizert2021b). Vertical gas diffusion within the firn pore network effectively stops at the lock-in depth L. At this depth, the gravitational isotopic enrichment of gases (such as δ¹⁵N-N₂) halts, and the Δage is fixed (Sowers and others, Reference Sowers, Bender, Raynaud and Korotkevich1992). The ice-equivalent lock-in depth (L _IE) gives the amount of ice (in meters) above L, and is given by

(1)\begin{equation}{L_{{\text{IE}}}} = \int \limits_0^L \rho \left( z \right)/{\rho _{{\text{ice}}}}dz\end{equation}

Over a wide range of climatic conditions, the ratio ${L_{{\text{IE}}}}/L$ is constant and ∼0.7 (Parrenin and others, Reference Parrenin2012). Presently at South Pole, this ratio is 0.712 (Vandecrux and others, Reference Vandecrux2023). The SPC accumulation rate is therefore given by

(2)\begin{equation}A = \,\frac{{{L_{{\text{IE}}}}}}{{{\Delta \text{age}}}} = \frac{L}{{0.712 \times {\Delta \text{age}}}}\end{equation}

For SPC, the Δage has been empirically reconstructed by combining volcanic synchronization of the ice ages, and methane synchronization of the gas ages to the WAIS Divide ice core that has an accurate layer-counted chronology and a small Δage (Buizert and others, Reference Buizert2015; Epifanio and others, Reference Epifanio2020; Sigl and others, Reference Sigl2016). In our reconstruction, we add 25 years to the published empirical Δage values to account for the diffusive age of the gases and thereby estimate the true ice age of lock-in (Buizert and others, Reference Buizert, Sowers and Blunier2013). Past L is reconstructed directly from the SPC δ¹⁵N data using the barometric equation (Sowers and others, Reference Sowers, Bender, Raynaud and Korotkevich1992), where we apply a cubic spline for smoothness and assume a constant 4 m convective zone thickness. We use standard error propagation techniques to estimate the uncertainty in A from the uncertainties in the empirical Δage (Epifanio and others, Reference Epifanio2020) and in L (set to 5 m). The resulting (1σ) uncertainty ranges from 5% to 15% of A. The A values reconstructed in this manner do not reflect the instantaneous values, but rather the values averaged over a time period corresponding the lifetime of the firn layer (ice equivalent, Δage). For the purpose of interpreting BNDs with respect to paleo-temperatures, this is actually the more meaningful parameter. Our accumulation record is shown in Fig. 5. Errors bars represent a root-mean-square (RMS) combination of the calculated gas-age/ice-age (Δage), and lock-in depth ( ${L_{{\text{IE}}}}$) uncertainties.

Figure 5. Firn-averaged accumulation (A) history (cm ice a⁻¹) for the South Pole site shown with combined uncertainty bands (gray) and discrete bubble number-density sampling depths (red). Published modern mean accumulation estimates are shown in ice-equivalent values.