The new Kr-86 excess ice core proxy for synoptic activity: West Antarctic storminess possibly linked to Intertropical Convergence Zone (ITCZ) movement through the last deglaciation

. Here we present a newly developed ice core gas-phase proxy that directly samples a component of the large-scale atmospheric circulation: synoptic-scale pressure variability. Surface pressure changes weakly disrupt gravitational isotopic settling in the ﬁrn layer, which is recorded in krypton-86 excess ( 86 Kr xs ). The 86 Kr xs may therefore reﬂect the time-averaged synoptic pressure variability over several years (site “storminess”), but it likely cannot record individual


Motivation and objectives
Proxy records from around the globe show strong evidence for past changes in Earth's atmospheric circulation and hydrological cycle that often far exceed those seen in the relatively short instrumental period.
For example, low-latitude records of riverine discharge captured in ocean sediments (Peterson et al., 2000) and isotopic composition of meteoric water captured in dripstone calcite (Cheng et al., 2016) suggest large variations in tropical hydrology and monsoon strength, commonly interpreted as meridional migrations of the Intertropical Convergence Zone or ITCZ (Chiang and Friedman, 2012;Schneider et al., 2014).Such ITCZ movement is seen in response to insolation changes linked to planetary orbit (Cruz et al., 2005) as well as in response to the abrupt millennial-scale Dansgaard-Oeschger (D-O) and Heinrich cycles of the North Atlantic (Kanner et al., 2012;Wang et al., 2001); the organizing principle is that the ITCZ follows the thermal Equator and therefore migrates towards the warmer (or warming) hemisphere (Broccoli et al., 2006;Chiang and Bitz, 2005).
As a second example, the intensity of the El Niño-Southern Oscillation (ENSO), the dominant mode of global interannual climate variability, has changed through time.A variety of proxy data suggest that ENSO activity in the 20th century was much stronger than in preceding centuries (Emile-Geay et al., 2015;Fowler et al., 2012;Gergis and Fowler, 2009;Thompson et al., 2013).The vast majority of data and model studies suggest weakened ENSO strength in the middle and early Holocene, likely in response to stronger orbitally driven NH summer insolation at that time (Braconnot et al., 2012;Cane, 2005;Clement et al., 2000;Driscoll et al., 2014;Koutavas et al., 2006;Liu et al., 2000Liu et al., , 2014;;Moy et al., 2002;Rein et al., 2005;Tudhope et al., 2001;Zheng et al., 2008); yet other studies suggest there may not be such a clear trend and simply more variability (Cobb et al., 2013).Intensification of ENSO (or perhaps a more El Niño-like mean state) may have occurred during the North Atlantic cold phases of the abrupt D-O and Heinrich cycles (Braconnot et al., 2012;Merkel et al., 2010;Stott et al., 2002;Timmermann et al., 2007).Overall, understanding past and future ENSO variability remains extremely challenging (Cai et al., 2015).
As a last example, the strength and meridional position of the Southern Hemisphere westerlies (SHWs) are thought to have changed in the past, which, via Southern Ocean winddriven upwelling, has potential implications for the global overturning circulation (Marshall and Speer, 2012) and for carbon storage in the abyssal ocean (Anderson et al., 2009;Russell et al., 2006;Toggweiler et al., 2006).The SHWs are thought to have been shifted equatorward (Kohfeld et al., 2013) during the Last Glacial Maximum (LGM), a shift on which climate models disagree (Rojas et al., 2009;Sime et al., 2013).During the abrupt D-O and Heinrich cycles, the SHWs move in parallel with the aforementioned migrations of the ITCZ in both data (Buizert et al., 2018;Marino et al., 2013;Markle et al., 2017) and models (Lee et al., 2011;Pedro et al., 2018;Rind et al., 2001).
As these examples clearly illustrate, evidence of past changes to the large-scale atmospheric circulation is widespread.However, proxy evidence of such past changes is typically indirect -for example via isotopes in precipitation, sea surface temperature, ocean frontal positions, windblown dust, or ocean upwelling -complicating their interpretation.
Here we present a newly developed noble-gas-based ice core proxy, Kr-86 excess ( 86 Kr xs ), that directly samples a component of the large-scale atmospheric circulation: synopticscale pressure variability.Owing to the firn air residence time of several years (Buizert et al., 2013) and the gradual bubble trapping process, each ice core sample contains a distribution of gas ages rather than a single age.Therefore, 86 Kr xs does not record the passing of individual weather systems, but rather the time-averaged intensity of synoptic-scale barometric variability.
Here we provide the first complete description of this new proxy.We validate and calibrate 86 Kr xs using late Holocene ice core samples from locations around Antarctica and Greenland that represent a wide range of pressure variability in the modern climate.We discuss the difficulties in using this proxy (analytical precision, surface melt, corrections for sample gas loss and thermal fractionation).Next, we use reanalysis data to better understand the drivers of surface pressure variability in Antarctica.Last, we present 86 Kr xs records from the Antarctic WAIS Divide ice core through the last deglaciation.

Gravitational disequilibrium and Kr-86 excess
The upper 50-100 m of the ice sheet accumulation zone consists of firn, the unconsolidated intermediate stage between snow and ice.An interconnected pore network exists within the firn, in which gas transport is dominated by molecular diffusion (Schwander et al., 1993).Diffusion in this stagnant air column results in gravitational enrichment in heavy gas isotopic ratios such as δ 15/14 N-N 2 , δ 40/36 Ar, and δ 86/82 Kr (Schwander, 1989;Sowers et al., 1992).In gravitational equilibrium, all these gases attain the same degree of isotopic enrichment per unit mass difference: with m the isotopic mass difference (1 × 10 −3 kg mol −1 ), g the gravitational acceleration, z the depth, R the gas constant, and T the temperature in Kelvin.
Besides molecular diffusion, firn air is mixed and transported via three other processes: downward advection with the sinking ice matrix, convective mixing (used in the firn air literature as an umbrella term to denote vigorous air exchange with the atmosphere via, e.g., wind pumping and seasonal convection), and dispersive mixing.These last three transport processes are all driven by large-scale air movement that does not distinguish between isotopologues, and we refer to them collectively as macroscopic air movement.Of particular interest for our proxy is dispersive mixing, which is driven by surface pressure variations.When a lowpressure (high-pressure) system moves into the site, firn air at all depth levels is forced upwards (downwards) to reach hydrostatic equilibrium with the atmosphere -a process called barometric pumping.One can think of the firn layer "breathing" in and out in response to a rising and falling barometer, respectively.Because firn has a finite dispersivity (Schwander et al., 1988), this air movement mixes the interstitial firn air.Note that an upward air movement also exists in the firn column relative to the overall downward advection of the ice, which is caused by the slow reduction of porosity with depth (Rommelaere et al., 1997).This upward airflow due to gradual pore closure (around 10 −9 to 10 −8 m s −1 ) is orders of magnitude smaller than the flows driven via barometric pumping (around 10 −6 m s −1 ) and therefore neglected here (Buizert and Severinghaus, 2016).
Any type of macroscopic air movement disturbs the gravitational settling, reducing isotopic enrichment below δ grav .Let δ 86 Kr, δ 40 Ar, and δ 15 N refer to deviations of 86 Kr/ 82 Kr, 40 Ar/ 36 Ar, and 29 N 2 / 28 N 2 , respectively, from their ratios in the well-mixed atmosphere.Gases that diffuse faster (such as N 2 ) will always be closer to gravitational equilibrium than gases that diffuse slower (such as Kr), and in the absence of thermal fractionation δ 86 Kr/4 < δ 40 Ar/4, δ 15 N < δ grav .The isotopic differences δ 86 Kr/4-δ 40 Ar/4 and δ 86 Kr/4-δ 15 N thus reflect the degree of gravitational disequilibrium.The magnitudes of the isotopic disequilibria scale in a predictable way following the molecular diffusion coefficients (Birner et al., 2018); because the diffusion coefficients of N 2 and Ar are very similar, their disequilibria are comparable in magnitude.We define Kr-86 excess using the Kr and Ar isotopic difference: where the "corr" subscript denotes a correction for gas loss (Appendix A1) and thermal fractionation (Appendix A2).
The rationale for including a normalization in the denominator is discussed below.An alternative Kr-86 excess definition is possible using δ 15 N instead of δ 40 Ar: Note that both definitions rely on having measurements of three isotope ratios (δ 86 Kr, δ 40 Ar, and δ 15 N), as the ther-mal correction requires δ 40 Ar and δ 15 N to be known.The 86 Kr xs40 definition is preferred because per unit mass difference δ 40 Ar is less sensitive to thermal fractionation than δ 15 N is (Grachev and Severinghaus, 2003a, b); this makes it more suitable for interpreting time series.Unless explicitly stated otherwise, we use 86 Kr xs40 as our definition of Kr-86 excess.The 86 Kr xs15 does provide a way to check the validity of 86 Kr xs40 time series, and indeed we find good correspondence between both definitions for the WDC deglacial time series (compare, for example, Figs.A3 and A4).Because the disequilibrium signal is small, we express 86 Kr xs in units of per meg (parts per million) of gravitational disequilibrium per per mil of gravitational enrichment.This unit (per meg ‰ −1 ) is mathematically identical to per mil, but we use it to emphasize the normalization in the denominator.
In the (theoretical) case of full gravitational equilibrium (and no gas loss or thermal fractionation), δ 86 Kr/4 = δ 40 Ar/4 = δ 15 N = δ grav , and therefore 86 Kr xs = 0. Any type of macroscopic mixing will cause δ 86 Kr/4 < δ 40 Ar/4 < δ 15 N < δ grav , and thus 86 Kr xs < 0. In this sense 86 Kr xs is a quantitative measure for the degree of gravitational disequilibrium in the firn layer (Birner et al., 2018;Buizert and Severinghaus, 2016).Kawamura et al. (2013) first describe this gravitational disequilibrium (or kinetic) fractionation effect at the Megadunes site (Severinghaus et al., 2010), where deep firn cracking leads to a 23 m thick convective zone.They suggest that the isotopic disequilibrium can be used to estimate past convective zone thickness.We show here that sites with small convective zones can nevertheless have very negative 86 Kr xs , and instead we suggest that the ice core 86 Kr xs is dominated by dispersive mixing driven by barometric pumping from timeaveraged synoptic-scale pressure variability.
The principle behind 86 Kr xs is illustrated with idealized firn model experiments in Fig. 1.In the absence of dispersive mixing (Fig. 1a, left panel), all isotope ratios approach δ grav and δ 86 Kr-δ 40 Ar is close to zero -but not exactly zero owing to downward air advection.Next, we replace a fraction f of the molecular diffusion with dispersive mixing.With dispersive mixing at f = 0.1 and f = 0.2 of total mixing (middle and right panels, respectively), isotopic enrichment is progressively reduced below δ grav (dashed line), making δ 86 Krδ 40 Ar (and consequently 86 Kr xs ) increasingly negative.
The ratio of macroscopic over diffusive transport is expressed via the dimensionless Péclet number, given here for advection and dispersion: where P e X is the Péclet number for gas X, w air the (downward) advective air velocity, L a characteristic length scale, D X the diffusion coefficient for gas X, and D disp the dispersion coefficient (Buizert and Severinghaus, 2016).In agreement with earlier studies (Birner et al., 2018;Kawamura et al., 2013), we find that δ 86 Kr-δ 40 Ar is maximized when molecular mixing and dispersive mixing are equal in magnitude (f = 0.5, Fig. 1b), corresponding to P e X ≈ 1.Note that 86 Kr xs responds more linearly to f than δ 86 Kr-δ 40 Ar does due to δ 40 Ar in the denominator of Eq. ( 2).In a last idealized experiment, we keep the fraction of dispersion fixed at f = 0.1 while we reduce the thickness of the firn column by increasing the site temperature (Fig. 1c).We find that δ 86 Kr-δ 40 Ar scales linearly with firn thickness, here represented by δ 40 Ar on the x axis.However, 86 Kr xs remains essentially constant due to the normalization by δ 40 Ar in the denominator of Eq. ( 2).The normalization step is thus necessary to enable meaningful comparison between different sites and time periods that all have different firn thicknesses.For this reason, the definition of 86 Kr xs used here has been updated from the original definition by Buizert and Severinghaus (2016).
Note that these highly idealized experiments assume dispersive mixing to be a fixed fraction of total transport throughout the firn column, equivalent to a constant Péclet number in the diffusive zone (a convective zone is absent in these simulations).In reality, the Péclet number varies greatly on all spatial scales.On the macroscopic scale (> 1 m), P e reflects the various transport regimes (Sowers et al., 1992), being highest in the convective and lock-in zones.On the microscopic scale (< 1 cm), hydraulic conductance scales as ∝ r 4 (with r the pore radius), whereas the diffusive conductance scales as ∝ r 2 .This means that the Darcy flow associated with barometric pumping will concentrate in the widest pores and pathways, leading to a range of effective Péclet numbers within a single sample of firn.At intermediate spatial scales of a few centimeters, firn density layering introduces strong heterogeneity in transport properties.It is unclear at present whether the competition between diffusive and non-diffusive transport, which occurs at the microscopic pore level, can be accurately represented in macroscopic firn air models via a linear parameterization as is the current practice.

Ice core sites
In this study we use ice samples from 11 ice cores drilled in Antarctica and 1 in Greenland.The Antarctic sites are the West Antarctic Ice Sheet (WAIS) Divide core (WDC06A, or WDC), Siple Dome (SDM), James Ross Island (JRI), Bruce Plateau (BRP), Law Dome DE08, Law Dome DE08-OH, Law Dome DSSW20K, Roosevelt Island Climate Evolution (RICE), Dome Fuji (DF), EPICA (European Project for Ice Coring in Antarctica) Dome C (EDC), and South Pole Ice Core (SPC14, or SP).Ice core locations in Antarctica are shown in Fig. 2a.In Greenland, we use samples from the Greenland Ice Sheet Project 2 (GISP2).
We shall refer to late Holocene data from these sites as the calibration dataset, analogous to a core-top dataset in the sediment coring literature.Site characteristics, coordinates, and number of samples included in the calibration dataset are given in Table 1.The DE08-OH site is a recent revisit of the Law Dome DE08 site.The DE08-OH core was measured at sub-annual resolution to understand centimeter-scale 86 Kr xs variations due to, for example, layering in firn density and bubble trapping (Appendix B).In addition to the calibration dataset, we present a record of Kr-86 excess going back to the LGM from WDC.

Ice sample analysis
We broadly follow analytical procedures described elsewhere (Bereiter et al., 2018a, b;Headly and Severinghaus, 2007;Severinghaus et al., 2003).In short, an 800 g ice sample, its edges trimmed with a band saw to expose fresh surfaces, is placed in a chilled vacuum flask that is then evacuated for 20 min using a turbomolecular pump.Air is extracted from the ice by melting the sample while stirring vigorously with a magnetic stir bar, led through a water trap, and cryogenically trapped in a dip tube immersed in liquid He.Next, the sample is split into two unequal fractions.The smaller fraction (about 2 % of total air) is analyzed for δ 15 N-N 2 , δ 18 O-O 2 , δO 2 /N 2 , and δAr/N 2 on a 3 kV Thermo Finnigan Delta V plus dual inlet IRMS (isotope ratio mass spectrometer).In the larger fraction, noble gases are isolated via hot gettering to remove reactive gases.The purified noble gases are then analyzed for δ 40/36 Ar, δ 40/38 Ar, δ 86/82 Kr, δ 86/84 Kr, δ 86/83 Kr, δKr/Ar, and δXe/Ar on a 10 kV Thermo Finnigan MAT253 dual-inlet IRMS.We reject one sample from RICE due to incomplete sample transfer and one sample from WDC due to problems with the water trap.Calibration is done for each measurement campaign by running samples of La Jolla pier air.
All calibration (core-top) data were measured using Method 2 as described by Bereiter et al. (2018a), with a longer equilibration time during the splitting step than used in that study to improve isotopic equilibration between the fractions.The exception is the DE08-OH site, where the ice sample (rather than the extracted gas sample) was split into two fractions -the advantage of this approach is that it does not require a gas splitting step that is time-consuming and may fractionate the isotopes; the downside is that the samples may have slightly different isotopic composition due to the stochastic nature of bubble trapping and the different gas loss histories of the ice pieces.
Measurements of the WDC downcore dataset were performed over five separate measurement campaigns that occurred in February-April 2014, February-April 2015, August 2015, August 2020, and August 2021.The first three campaigns are described by Bereiter et al. (2018b), in which the 86 Kr xs data are a by-product of measuring δKr/N 2 for reconstructing global mean ocean temperature.Campaigns 1 and 2 are in good agreement, whereas campaign 3 appears offset from the other two by an amount that exceeds the analytical precision (offset around 35 per meg ‰ −1 ).To validate the main features in the record, we performed two additional campaigns (4 and 5), in which all the gas extracted from each ice sample was quantitatively gettered and only analyzed for Ar and Kr isotopic composition.The downcore record, as well as the five analytical campaigns, is discussed in detail in Sect.5.1.Data from the bubble-clathrate transition zone (here 1000 to 1600 m depth, or ∼ 4 ka to 7 ka) are excluded owing to the potential for artifacts; the depth range of the bubble-clathrate transition zone is based on observed positive anomalies in δO 2 /N 2 in WDC ice.
All samples were analyzed at the Scripps Institution of Oceanography, USA, with the exception of the EDC samples, which were analyzed at University of Bern, Switzerland    (Baggenstos et al., 2019).Some of the EDC samples analyzed had clear evidence of drill liquid contamination, which acts to artifactually lower 86 Kr xs via isobaric interference on mass 82; the late Holocene data used here were not flagged for drill liquid contamination (Baggenstos et al., 2019).The 2σ analytical precision of the δ 15 N, δ 40 Ar, and δ 86 Kr measurements is around 3, 5, and 26 per meg, respectively, based on the reproducibility of La Jolla air measurements.Via standard error propagation, this results in a ∼ 22 per meg ‰ −1 (2σ ) analytical uncertainty for both 86 Kr xs40 and 86 Kr xs15 at a site like WDC where δ 40 Ar ≈ 1.2 ‰.We have no true (same-depth) replicates to assess the reproducibility of 86 Kr xs measurements experimentally.The measured isotope ratios are corrected for gas loss ( 40 GL ) and thermal fractionation ( 86 TF , 40 TF , 15 TF ) before interpretation; details on these corrections are given in Appendix A. For the coretop calibration study, the average magnitude of the gas loss and thermal fractionation corrections is +14 and −15 per meg ‰ −1 in 86 Kr xs , respectively.Note that these two corrections both involve the δ 40 Ar isotopic ratio, and therefore they are not independent from each other and not additivein other words, the total correction is not simply the sum of the two individual corrections.
Our study includes two ice cores from the Antarctic Peninsula: BRP (two ice samples) and JRI (five ice samples).Measured δXe/N 2 ratios (and to a lesser extent the δKr/N 2 ratios) in all samples from both locations are significantly elevated above the expected gravitational enrichment signal (Fig. A1a), which is clear evidence for the presence of refrozen meltwater in these samples (Orsi et al., 2015).Like xenon, krypton is highly soluble in (melt)water, and therefore 86 Kr xs cannot be reliably measured in these samples; we reject all samples from the BRP and JRI sites.It is notable that all samples from both sites show evidence of refrozen meltwater, given that the high-accumulation BRP core is nearly entirely free of visible melt layers and that we carefully selected samples without visible melt features at JRI. Visible ice lenses form only when meltwater pools and refreezes on top of low-permeability layers such as wind crusts; our ob-servations suggest meltwater can also refreeze throughout the firn in a way that cannot be detected visually.

Calibrating Kr-86 excess
The 86 Kr xs proxy for synoptic activity was first proposed on theoretical grounds by Buizert and Severinghaus (2016) -here we provide the first experimental validation of this proxy using a core-top calibration of 86 Kr xs using late Holocene ice core samples from nine locations around Antarctica and one in Greenland that represent a wide range of pressure variability in the modern climate (here: 1979-2017 CE).

Spatial variation in synoptic-scale pressure variability
Kr-86 excess is sensitive to air movement (both upward and downward), which in turn is controlled by the magnitude of relative air pressure change.Let p i be a time series of (synoptic-scale) site surface pressure with N data points, time resolution t, and mean value p.The time series can span a month, year, or multi-year period, with p potentially different for each month or year.We define the parameter as which we express here in convenient units of percent per day (% d −1 ).reflects the intensity of barometric pumping in the firn column.Note that t should be larger than ∼ 1 h, which is the timescale for the entire firn column to equilibrate with the surface pressure (Buizert and Severinghaus, 2016), and smaller than about a day in order to adequately resolve synoptic-scale pressure events.Here we use ERA-Interim reanalysis data from 1979-2017 with t = 6 h (Dee et al., 2011), from which we calculate monthly and annual values using Eq. ( 5).A map of annual mean across Antarctica is given in Fig. 2a.At all sites considered, has a strong seasonal cycle with pressure variability and storminess being strongest in the local winter season (Fig. 2c).Interannual differences in are greatest along the Siple coast and coastal West Antarctica (Fig. 2b), mainly reflecting the influence of central Pacific (ENSO, Pacific Decadal Oscillation -PDO) climate variability (Sect.4).

Kr-86 excess proxy calibration
Present-day Antarctica has a wide range of (Fig. 2a), which allows us to validate and calibrate 86 Kr xs .In Fig. 3a we plot the site mean 86 Kr xs40 (with ±1σ error bars) as a function of (averaged over full 1979-2017 period).We find a Pearson correlation coefficient of r = −0.94when using site mean 86 Kr xs40 and r = −0.83when using the 86 Kr xs40 of individual samples (p < 0.01).Note that in this particular case the site mean 86 Kr xs40 and 86 Kr xs15 are identical (because by design, after thermal correction δ 15 N = δ 40 Ar/4); the error bars are different, though.
The 86 Kr xs data have been corrected for gas loss (Appendix A1) and thermal fractionation (Appendix A2), with the gas loss correction being the more uncertain component.Figure 3b shows the correlations of the calibration curve as a function of the gas loss scaling parameter ε 40 .We find a good correlation over a wide range of ε 40 values, proving our calibration is not dependent on the choice of ε 40 .When using uncorrected 86 Kr xs40 data the site mean correlation is r = −0.71;when applied individually, both the gas loss and thermal correction each improve the correlation to r = −0.77and r = −0.79,respectively (Fig. A3, all p < 0.05).Based on these tests we conclude that the observed relationship is not an artifact of the applied corrections.The applied corrections improve the correlation, which increases confidence in the method.The calibration results for 86 Kr xs15 are shown in Fig. A4.
Notably, there is a large spread in 86 Kr xs across samples from any single site, particularly at the high-sites of SDM and RICE (Fig. 3a, note the ±1σ error bars).This spread is larger than the measurement uncertainty, and we believe this variance reflects a signal that is truly present in the ice.The Siple coast and Roosevelt Island experience the largest interannual differences in Antarctica (Fig. 2b), and it is therefore likely that our coarse sampling is aliasing the true 86 Kr xs signal.The variance in 86 Kr xs may also contain climate information.
Both theoretical considerations and observations thus suggest 86 Kr xs is a proxy for time-averaged barometric surface pressure variability at the site, and in the remainder of this paper we will interpret it as such.

Discussion of the Kr-86 excess proxy
Our interpretation of 86 Kr xs as a proxy for time-averaged pressure variability is somewhat complicated by the possibility of deep convective zones, which have the same 86 Kr xs signature as barometric pumping.This was discovered at the Megadunes (MD) site, central East Antarctica; at this zeroaccumulation site deep cracks form in the firn layer that facilitate a 23 m deep convection zone (Severinghaus et al., 2010).In fact, this observation led earlier work to suggest that noble gas gravitational disequilibrium may be used as a proxy for convective zone thickness (Kawamura et al., 2013) rather than synoptic-scale pressure variability as suggested here.Although megadunes and zero-accumulation zones are ubiquitous and cover 20 % of the Antarctic Plateau (Fahnestock et al., 2000), ice cores are seldom drilled in these areas and it is safe to assume that they never formed at sites like WAIS Divide that had relatively high accumulation rates even during the last glacial period.Performing the corrections for thermal and size-dependent fractionation is challenging at MD, and we suggest that the MD 86   Triangles denote the ε 40 estimate from the Byrd, Siple, and GISP2 ice cores (Fig. A2; Kobashi et al., 2008a;Severinghaus et al., 2003).
−55 per meg ‰ −1 ; even at the larger limit, this is still smaller in magnitude than 86 Kr xs anomalies at several modern-day sites with small convective zones (such as SDM, RICE, and the Law Dome sites), suggesting barometric pumping is capable of producing larger 86 Kr xs signals than even the most extreme observed case of convective surface mixing.Having 86 Kr xs measured in MD ice core (rather than firn air) samples would be valuable for a more meaningful comparison to the ice core sample measurements presented here.Windy sites can have substantial convective zones of ∼ 14 m (Kawamura et al., 2006), and future studies of 86 Kr xs at such sites would be valuable.
Currently, 1-D and 2-D firn air transport model simulations underestimate the magnitude of the 86 Kr xs signal compared to measurements in mature ice samples (Birner et al., 2018), complicating scientific understanding of the proxy.In these models, the effective molecular diffusivity of each gas is scaled linearly to its free air diffusivity.The ratio of krypton to argon free air diffusivity is 0.78.This ratio, which directly sets the magnitude of the simulated 86 Kr xs , may actually be smaller than 0.78 in real firn, as krypton is more readily adsorbed onto firn surfaces, retarding its movement (similar to gases moving through a gas chromatography column).This may be one explanation for why models simulate too little 86 Kr xs .
Another likely explanation for the model-data mismatch is that certain critical sub-grid processes (such as the aforementioned pore size dependence of the Péclet number) are not adequately represented in these models.Barometric pumping may further actively shape the pore network through the movement of water vapor, thereby keeping certain preferred pathways connected and open below the density at which percolation theory would predict their closure (Schaller et al., 2017).The fate of a pore restriction is determined by the balance between the hydrostatic pressure (that acts to close it) and vapor movement away from its convex surfaces (that acts to keep it open); we speculate that barometric Darcy airflow keeps high-flow channels connected longer by eroding convex surfaces.This enhances the complexity (and therefore dispersivity) of the deep firn pore network and possibly creates a nonlinear 86 Kr xs response to barometric pumping.The hypothesized channel formation in deep firn is driven by a positive feedback on flow volume and somewhat reminiscent of erosion-driven stream network formation in fluvial geomorphology.
Firn models predict that, after correcting for thermal fractionation, the deviation from gravitational equilibrium for the elemental ratios (such as δKr/Ar) should be proportional to that deviation in isotopic ratios.However, the observations suggest that the former is usually smaller than would be expected from the latter.As before, adsorption of Kr onto firn grain surfaces may contribute to the observed discrepancy, and laboratory tests of this process are called for.Further, the impacts of gas loss are greater on elemental ratios than on the isotopic ratios, which may also contribute.Including measurements of xenon isotopes and elemental ratios in future measurement campaigns may be able to provide additional constraints to better understand this discrepancy.
Measurements on firn air samples, where available, suggest a smaller 86 Kr xs anomaly in firn air than found in ice core samples from the same site.We attribute this in part to a seasonal bias that is introduced by the fact that firn air sam-pling always takes place during the summer months, whereas the synoptic variability that drives the Kr-86 excess anomalies is largest during the winter (Fig. 2c); consequently, firn air observations are biased towards weaker 86 Kr xs .Further, in the deep firn where 86 Kr xs anomalies are largest, firn air pumping may not yield a representative air sample, but rather be biased towards the well-connected porosity at the expense of poorly connected cul-de-sac-like pore clusters.Since barometric pumping ventilates this well-connected pore space with low-86 Kr xs air from shallower depths, the firn air sampling may not capture a representative 86 Kr xs value of the full firn air content.These explanations are all somewhat speculative, and a definitive understanding of the firn-ice differences is lacking at this stage.
Gas loss and thermal corrections are critical to the interpretation of 86 Kr xs .The thermal correction is applied to account for thermal gradients in the firn ( T , defined here as the temperature at the top minus the temperature at the base of the firn), which are chiefly caused by geothermal heat or surface temperature changes at the site.At low-accumulation sites geothermal heating leads to T < 0. We use 15 N excess (δ 15 N-δ 40 Ar/4) to estimate the thermal gradient in the firn (Appendix A2).Because nitrogen and argon have similar diffusivities but different thermal diffusion coefficients, δ 15 Nδ 40 Ar is relatively insensitive to barometric pumping yet sensitive to thermal fractionation, allowing estimating T .
Besides the actual thermal gradients in the firn, the isotopic composition may also be impacted by seasonal rectifier effects.If the firn air transport properties differ between the seasons (for example due to thermal contraction cracks, convective instabilities, or seasonality in wind pumping), this can result in a thermal fractionation of isotopic ratios in the absence of a thermal gradient T in mean annual temperature (Morgan et al., 2022).
For the WDC, DSS, and GISP2 sites we obtain T values close to zero as expected for these high-accumulation sites; for the SP, SDM, RICE, and DF sites we find T ranging from −0.76 to −1.18 • C, in agreement with the effect of geothermal heat.The high-accumulation DE08 and DE08-OH sites both have an unexpectedly large T of −1.6 • C; the good agreement between the sites suggests that it is likely a real signal, yet we can rule out geothermal heat as the cause.This may suggest that the Law Dome DE08 site is subject to a seasonal rectifier effect or a recent climatic cooling.Last, the EDC site shows an unexpected positive T = +1.6 ± 1.89 • C. Three possible explanations are the following (1) the aforementioned drill liquid contamination for this core (Baggenstos et al., 2019); (2) a summertimebiased seasonal rectifier; or (3) an overcorrection of δ 40 Ar for gas loss, which could occur, for example, if natural gas loss and post-coring fugitive gas loss fractionate δ 40 Ar differently and EDC samples were impacted mostly by the former type (our correction is mostly based on measurements of the latter type).
For the Law Dome DE08-OH site we observe large (5fold) sub-annual variations in 86 Kr xs (Fig. B1).The magnitude of the 86 Kr xs layering is truly remarkable.The isotopic enrichment of each gas (δ 15 N, δ 40 Ar, δ 86 Kr) can be converted to an effective diffusive column height (DCH).For the samples with the smallest (greatest) 86 Kr xs magnitude, this DCH is around 1 m (6 m) shorter for δ 86 Kr than it is for δ 15 N.The firn air transport physics that may explain such phenomena are beyond our current scientific understanding.The sub-annual variations may be related to the seasonal cycle in storminess, though that seems improbable to us at present as the gas age distribution at the depth of bubble closure has a width of several years (Schwander et al., 1993).Another reason may be seasonal layering in firn properties -such as density, grain size, and pore connectivity -that control the degree of disorder and dispersive mixing occurring in the firn and lead to a staggered firn trapping and seasonal variations in age (Etheridge et al., 1992;Rhodes et al., 2016).The sample air content estimated from the IRMS inlet pressure is similar for all measurements, making it unlikely that the variations in 86 Kr xs are caused by remnant open porosity in lower-density layers.In any case it is remarkable that such large variations in gas composition can arise and persist on such small length scales, given the relatively large diffusive, dispersive, and advective transport length scales of the system.More work is needed to establish the origin of the sub-annual variations in ice core 86 Kr xs .At all other sites analyzed here, the sample length exceeds the annual layer thickness; this will remove some, but not all, of the effects of the sub-annual variations.
Another puzzling observation is the positive 86 Kr xs at the Dome Fuji (DF) site; theoretical considerations suggest it should always be negative.In part this may be due to an overcorrection of δ 40 Ar for gas loss, which would act to bias 86 Kr xs in the positive direction.This correction is largest at DF owing to the very negative δO 2 /N 2 and δAr/N 2 (Fig. A1); while we base our correction on published work, it is conceivable that we overestimate the true correction (Appendix A1).In particular, our gas loss correction is based on observations of artifactual post-coring gas loss, which may fractionate δ 40 Ar differently than natural fugitive gas loss during bubble close-off.Omitting the gas loss correction indeed makes 86 Kr xs at DF negative (Fig. A3c and d).Another hypothesis is that the positive 86 Kr xs signal is an artifact of the seasonal rectifier that Morgan et al. (2022) identify at DF.In this work we assume a linear approach in which the effect of the rectifier can be described by a single T value that is the same for isotopic pairs.In reality, there may be nonlinear interactions between thermal fractionation and firn advection that impact the isotopic values of the various gases in a more complex way than captured in our approach.
The 86 Kr xs is also correlated with other site characteristics besides .For site elevation we find r = 0.96 (0.84) and for mean annual temperature r = −0.87(−0.76); the number in parentheses gives the correlation when using all the https://doi.org/10.5194/cp-19-579-2023Clim. Past, 19, 579-606, 2023 individual samples rather than site mean 86 Kr xs .The listed correlations all have p < 0.01.For site accumulation we do not find a statistically significant correlation at the 90 % confidence level.The correlations with elevation and temperature are comparable to those we find for ; this is no surprise given that elevation, , and T are all strongly correlated with one another, mainly because elevation directly impacts both T (via the lapse rate) and (via its topographic influence on the position of storm tracks).To our knowledge there are no mechanisms through which either elevation or annual mean temperature could drive kinetic isotopic fractionation in the firn layer.Perhaps other unexamined site characteristics (such as the degree of density layering or the magnitude of the annual temperature cycle) could also provide good correlations, suggesting additional hidden controls on 86 Kr xs .The data needed to assess such hidden controls are not available for most sites.
The calibration of the 86 Kr xs proxy is based on spatial regression.In applying the proxy relationship to temporal records, we make the implicit assumption that proxy behavior in the temporal and spatial dimensions is at least qualitatively similar.This assumption may prove incorrect.In particular, changes in insolation are known to impact firn microstructure and bubble close-off characteristics, which in turn impacts gas records of δO 2 /N 2 and total air content (Bender, 2002;Raynaud et al., 2007).Since 86 Kr xs is linked to the dispersivity of deep firn, it seems probable that insolation also has a direct impact on 86 Kr xs via the firn microstructure.We will revisit this issue in our interpretation of the WDC 86 Kr xs record (Sect.5).Overall, we anticipate 86 Kr xs to be a qualitative proxy for synoptic variability yet want to caution against quantitative interpretation based on the spatial regression slope.
The observations presented in this section clearly highlight the fundamental shortcomings of our current understanding of firn air transport, hinting at the existence of complex interactions, presumably at the pore scale, that are not being represented.Percolation theory finds that near the critical point (presumably the lock-in depth) a network becomes fractal in its nature; we suggest that this fractal nature of the pore network likely contributes to nonlinear pore-scale interactions that give rise to the 86 Kr xs observations in ice.While the observed correlation in Fig. 3a is highly encouraging, further work is critical to understand this proxy.Examples of such future studies are the following: (1) additional highresolution records that can resolve the true variations that exist in a single ice core, similar to the DE08-OH record; (2) 3-D firn air transport model studies; (3) improvements to the gas loss correction; (4) additional coring sites to extend the spatial calibration and further confirm the validity of the proxy; (5) adding xenon isotopic constraints ( 136 Xe excess) as an additional marker of isotopic disequilibrium; (6) numerical simulations of pore-scale air transport in large-scale firn networks; (7) experimental studies of dispersion and noble gas adsorption in firn samples; (8) percolation theory ap-proaches to study the fractal nature of the pore network of the lock-in zone; and (9) replication of the WDC deglacial 86 Kr xs record in nearby ice cores such as RICE.

Present-day controls on Kr-86 excess in Antarctica
In this section we investigate the large-scale patterns of climate variability in the Southern Hemisphere that could affect and therefore 86 Kr xs over Antarctica.We begin by investigating the patterns in the wind field that are associated with changes in at ice core sites, before examining how more canonical patterns of Southern Hemisphere climate variability, such as the southern annular mode (SAM), might affect over the whole of Antarctica.
We use ERA-Interim reanalysis data for the 1979-2017 period (Dee et al., 2011) to evaluate the present-day controls on synoptic-scale pressure variability in Antarctica.Kr-86 excess in an ice core sample averages over several years of pressure variability, and therefore we focus on annual mean correlation in our analysis.The annual mean is calculated from the 6-hourly reanalysis data using Eq. ( 5).Note that we let the year run from April to March to avoid dividing single El Niño and La Niña events across multiple years.
At all Antarctic sites investigated, a similar pattern exists; four representative locations are shown in Fig. 4, where we regress the zonal wind in the lower (850 hPa, color shading) and upper troposphere (200 hPa, contours) onto our surface pressure variability parameter .We find that synoptic pressure variability at these sites is linked to zonal winds along the southern margin of the eddy-driven subpolar jet (SPJ), which extends from the surface to the upper troposphere (Nakamura and Shimpo, 2004;Trenberth, 1991).Sites near the ice sheet margin (Fig. 4a, b and d) are most sensitive to the SPJ edge in their sector of Antarctica, whereas interior sites (Fig. 4c) appear to be sensitive to the overall strength and position of the SPJ.Note that strengthening, broadening, or southward shifting of the SPJ can all in principle enhance site .
Pressure variability at WDC is furthermore correlated with the strength of the Pacific subtropical jet (STJ) aloft (solid contour lines centered around 30 • S in the Pacific in Fig. 4a), forming an upper troposphere wind pattern that resembles the wintertime South Pacific split jet (Bals-Elsholz et al., 2001;Nakamura and Shimpo, 2004); this agrees with the finding that a strengthening of the split jet enhances storminess over West Antarctica (Chiang et al., 2014).
Next, we investigate how the well-known patterns of largescale atmospheric variability, such as SAM and ENSO, impact pressure variability in Antarctica. Figure 5 shows the correlation of with the three leading modes of SH extratropical atmospheric variability; the correlation with various indices and modes for individual ice core locations is given in Table 2. Most teleconnection patterns have a specific season during which they are strongest; here we do not differentiate between seasons because 86 Kr xs in ice core samples averages over all seasons.
Globally, annual mean is highest over the Southern Ocean (Fig. 5a), a region of enhanced baroclinicity associated with the eddy-driven SPJ (Nakamura and Shimpo, 2004).The green line denotes the latitude of maximum , roughly corresponding to the latitude with the highest storm track density (57.8 • S on average).
The dominant mode of atmospheric variability in the SH extratropics is the southern annular mode, representing the vacillation of atmospheric mass between the middle and high latitudes (Thompson and Wallace, 2000).Figure 5b shows 500 hPa geopotential height (Z 500 ) anomalies associated with the SAM as contours, with the color shading giving the correlation between and the SAM index.During the positive SAM phase (negative Z 500 over Antarctica) we find that the storm tracks and maximum synoptic activity are displaced towards Antarctica (positive correlation poleward of the green line in Fig. 5b).This is associated with a strengthening and poleward displacement of the SH westerly winds that occurs during a positive SAM phase.More locally, on the Antarctic Peninsula is positively correlated with the SAM index (Table 2); at the other sites is not meaningfully impacted.This suggests that the variations associated with the SAM (as commonly defined) do not extend far enough poleward to meaningfully impact Antarctica with the exception of the peninsula.Enhanced synoptic variability on the peninsula during positive SAM phases is consistent with observations of enhanced snowfall at those times (Thomas et al., 2008).
The second mode of SH extratropical variability is the Pacific-South American Mode 1 (PSA1), which reflects a Rossby wave response to sea surface temperature (SST) anomalies over the central and eastern equatorial Pacific (Mo and Paegle, 2001) and is therefore closely linked to ENSO on interannual timescales (we find a correlation of r = 0.77 between the annual mean PSA1 and Niño 3.  Table 2. Pearson correlation between at the ice coring sites and large-scale atmospheric circulation.Correlations are calculated using annual mean data (all months, April-March).We only list the statistically significant correlations (p < 0.1).The Niño 3.4 is calculated over 5 • S-5 and Ross Sea sectors (Fig. 5c and Table 2).While all the correlations listed are statistically significant, they explain only a fraction of the total variability.Next, we consider anomalies in sea ice area and extent (Parkinson and Cavalieri, 2012).We focus on the Ross and Amundsen-Bellingshausen seas where impacts on the WAIS Divide may be expected.At the 90 % confidence level we do not find significant correlations with sea ice area or extent at most core locations (Table 2).Correlations with sea ice extent are (even) weaker than those for sea ice area and consequently not shown.We performed a lead-lag study of the correlations between and sea ice area as well as extent in the various sectors and find that in all cases maximum correlations occur for the sea ice changes lagging 0 to 4 months behind ; we interpret this to mean that the sea ice is responding to changes in atmospheric circulation rather than driving them.
Overall, we find that synoptic activity at WAIS Divide, the site of most interest here, is controlled by the position and/or strength of the storm tracks at the southern edge of the SPJ in the Pacific sector of the Southern Ocean (Ross, Amundsen, and Bellingshausen seas), with little sensitivity to the SPJ behavior in the other sectors.Owing to its remote southern location, WDC is only weakly impacted by the commonly defined large-scale modes of atmospheric variability.Most notably, WDC has a modest influence from the tropical Pacific climate, as shown by a correlation around r ≈ 0.3 with the PSA1, Niño 3.4, and PDO indices (Table 2).We further find statistically significant correlations (up to r = 0.44) between WDC and SST in broad regions of the central and eastern tropical Pacific (not shown).We suggest that ENSO weakly impacts storminess at WDC (around 10 % of variance explained) via its impact on the SPJ in the South Pacific.

Barometric variability in West Antarctica during the last deglaciation
5.1 The 0-24 ka WAIS Divide Kr-86 excess record The WAIS Divide downcore 86 Kr xs dataset we present here was produced during five separate measurement campaigns that occurred in February-April 2014, February-April 2015, August 2015, August 2020, and August 2021.Campaigns 1-3 were reported previously (Bereiter et al., 2018a, b), and campaigns 4 and 5 were meant to resolve conflicts between the 86 Kr xs datasets from these earlier campaigns.Three slightly different measurement approaches were used.Campaign 1 uses Method 1 from Bereiter et al. (2018a), in which the air sample splitting is done in a water bath for over 12 h to equilibrate the sample.Campaigns 2 and 3 use Method 2 from Bereiter et al. (2018a), in which a bellows is used to split the air samples for over 4 to 6 h.Campaigns 4 and 5 do not involve splitting of the air sample and only analyzed the Kr and Ar isotopic ratios.During campaign 4 a glass bead from the water trap got stuck in the tubing, restricting the flow and likely resulting in incomplete air extraction from the meltwater.Figure 6 compares 86 Kr xs40 (Fig. 6a) and 86 Kr xs15 (Fig. 6b) from the five campaigns.Campaign 1 is the only campaign that spans the full age range of the record, making it the most valuable of the three campaigns.Campaigns 2 and 3 are mostly restricted to the Pleistocene and Holocene periods, respectively, with little overlap between them.Campaigns 4 and 5 aimed to reproduce some of the most salient features in the earlier three.For ease of visual comparison, a Gaussian smoothing spline to the combined data from the two most extensive campaigns (1 and 2) is added (details below).
Visual inspection suggests that campaigns 1 and 2 agree well for both 86 Kr xs definitions.Campaign 3 has more scatter and is visibly offset (more negative) for both 86 Kr xs definitions.Campaign 4 has several data points that agree well with the spline, yet also several fliers with very negative 86 Kr xs40 .Campaign 5 appears to be overall more negative in 86 Kr xs40 .No true replicate samples were analyzed between the campaigns, in part because the large sample size requirement precludes this.This precludes a direct assessment of campaign offsets.Instead we rely on linear interpolation.First, we linear interpolate the 86 Kr xs values of campaign 1 onto the sample depths of campaign 2, and vice versa, allowing us to estimate the offsets.In this way, during their period of overlap we estimate offsets of 6 and 13 per meg ‰ −1 between campaigns 1 and 2 for 86 Kr xs40 and 86 Kr xs15 , respectively.This is within the analytical precision (22 per meg), suggesting these two campaigns are in good agreement.Data from campaign 2 appear to have more scatter, possibly reflecting the shorter equilibration time during sample splitting.
We combine data from the first two campaigns and evaluate their offset to data from the other three campaigns again using the linear interpolation method.For campaigns 3, 4, and 5 we find an offset of −32, −22, and −23 per meg ‰ −1 in 86 Kr xs40 , respectively.For campaign 3 the offset is −34 per meg ‰ −1 in 86 Kr xs15 .It is remarkable that all three later campaigns are more negative in 86 Kr xs than the first two.Campaign 3 shows the greatest offset (greater than analytical precision) and has more scatter in both 86 Kr xs (Fig. 6) and 15 N excess, possibly because for this campaign less care was taken that the IRMS conditions were stable.The offset of campaign 4 may be attributed to the incomplete sample transfer due to the bead stuck in the line; note that for this campaign the offset is caused by two very negative data points.The offset in campaign 5 is hard to explain.The more negative 86 Kr xs of campaigns 4 and 5 may reflect sample storage effects, as these were measured 5-6 years after campaign 1 and 2. However, this would not explain the negative values of campaign 3. The good 86 Kr xs agreement between DE08 and DE08-OH, drilled 32 years apart, would also argue against large storage effects.For campaigns 4 and 5 only Ar and Kr isotope ratios were measured, so we lack typical tracers of gas loss (δO 2 /N 2 and δAr/N 2 ) that can be analyzed.
In the remainder of this paper we will interpret the combined data from campaigns 1 and 2, but with the caveat that there is a persistent offset with later campaigns.However, the features we interpret are corroborated by the later campaigns if one takes the offset into account.To aid interpretation of the data, we apply a Gaussian smoothing spline with a smoothing filter width that varies depending on the data density (from 250-year width in the deglaciation itself for which the data density is high to 1750 years in the Holocene and LGM for which data density is low).certainty in the smoothing spline we use a Monte Carlo approach that considers uncertainty in following: (1) the gas loss correction by randomly sampling ε 40 in the range of 0 to −0.016 as well as by randomly adding an offset in the range of −1 ‰ to +1 ‰ to the gas loss indicator (δO 2 /N 2 -δAr/N 2 ), (2) the thermal correction by randomly scaling the thermal scenario (Fig. A5) by a factor ranging from 0 to 2, and (3) analytical errors by adding random errors to individual data points drawn from a normal distribution with a 2σ width of 22 per meg.The ±1σ uncertainty range and mean value are shown as the gray envelope and center line in Fig. 6.We the following observations to be robust.
-The Holocene shows a trend towards increasingly negative 86 Kr xs , suggesting a gradual increase in synoptic activity toward the present.Minimum synoptic activity in West Antarctica occurs during the early Holocene around 10 ka; the Monte Carlo study suggests 86 Kr xs40 in the early Holocene (8 ka-10 ka) is 30.5 ± 18 per meg ‰ −1 (±2σ ) below the late Holocene value (last 2 ka).Using the slope of our core-top calibration (Fig. 3), we estimate that early Holocene WDC synoptic activity is ∼ 17 % weaker than it is today.This change is comparable to the 2σ magnitude of interannual variations in annual mean at the site today (or about half the peak-to-peak variations thereof).This Holocene trend is seen in the data from campaigns 1, 3, and 4. Campaign 5 does not suggest a trend but has only one late Holocene data point, making it less robust.The trend in campaign 3 is less robust due to the greater scatter in the data.
-The most pronounced change occurs at the Younger Dryas (YD)-Holocene transition, where 86 Kr xs becomes more positive (by 30.1 ± 16 per meg ‰ −1 , comparing YD and early Holocene), implying a decrease in synoptic activity.This transition is observed in campaigns 1, 2, 4, and 5 that cover this time period (the third campaign does not cover it) and represents a ∼ 17 % drop in synoptic activity ( ).
-During the Last Glacial Maximum (LGM), WDC synoptic activity was perhaps slightly weaker than at present, but not significantly so ( 86 Kr xs40 more positive by 11 ± 13 per meg ‰ −1 ).The West Antarctic Ice Sheet elevation was likely higher during the LGM, and a 300 m elevation increase would by itself increase 86 Kr xs40 by 10 per meg ‰ −1 , all else being equal (Appendix A3); this is within the analytical error of our observations.This feature is seen in campaign 1 and not covered by the other campaigns.
-The deglaciation itself has enhanced synoptic activity, in particular during the two North Atlantic cold stages Heinrich Stadial 1 (HS1) and the YD as highlighted with yellow bars in Figs. 6 and 7. Synoptic activity during these periods is enhanced relative to the adjacent LGM and early Holocene yet comparable to today.This feature is seen in campaigns 1 and 2 as well as in 4 and 5 for the transition into the Holocene.
Below we will interpret the deglacial WDC 86 Kr xs record in terms of time-averaged barometric variability.Before doing so we want to emphasize that firn processes may have also been imprinted onto the record, in particular on orbital timescales on which firn microstructure responds to local (summer) insolation intensity (Bender, 2002).High summer insolation results in more depleted δO 2 /N 2 and reduced air content, likely via stronger layering and a delayed pore closeoff process (Fujita et al., 2009).
Local summer solstice insolation in Antarctica increases through the Holocene, with the highest values in the late Holocene.This may impact 86 Kr xs , although it is not a priori clear what the sign of this relationship would be.The sense of the Holocene temporal trends is that a more negative 86 Kr xs coincides with more negative δO 2 /N 2 .Note that this is opposite to the trends seen in the spatial calibration, for which sites with the most negative δO 2 /N 2 (DF, SP, EDC) have the most positive 86 Kr xs .For now, the impact of local insolation on 86 Kr xs via firn microstructure remains unknown, which is an important caveat in interpreting the orbital-scale changes in WDC 86 Kr xs .The abrupt 86 Kr xs increase at the Holocene  (Grootes et al., 1993).(b) Hemispheric temperature difference (McGee et al., 2014) based on global proxy compilations for the Holocene (Marcott et al., 2013) and last deglaciation (Shakun et al., 2012).(c) Speleothem calcite δ 18 O from Hulu and Dongge caves, China, as a proxy for East Asian summer monsoon strength (Dykoski et al., 2005;Wang et al., 2001).Superimposed is summer solstice (21 June) insolation at 30 • N. (d) Speleothem calcite δ 18 O from Botuvera cave, southern Brazil, as a proxy for South American summer monsoon strength (Cruz et al., 2005;Wang et al., 2007).(e) Kr-86 excess record from WAIS Divide (this study), corrected for gas loss and thermal fractionation (Appendix A).The center line and shaded envelope show the mean and ±1σ uncertainty interval of a 10 000-iteration Monte Carlo smoothing exercise (see text).The dotted red line equals the center line with a correction for elevation change applied (Appendix A) using a simulated elevation history (Golledge et al., 2014).(f) Number of El Niño events per century from laminations in sediments from Laguna Pallcacocha, Ecuador (Moy et al., 2002).(g) Th-normalized opal flux in the Pacific Antarctic zone (south of the polar front) from cores NBP9802-6PC1 (turquoise; 169.98 • W, 61.88 • S) and PS75/072-4 (blue; 151.22 • W, 57.56 • S), reflecting local productivity and (wind-driven) upwelling (Chase et al., 2003;Studer et al., 2015).All water isotope data in this figure are on the Vienna Standard Mean Ocean Water (V-SMOW) scale.Arrows show the direction of increased monsoon strength and synoptic activity.
https://doi.org/10.5194/cp-19-579-2023Clim. Past, 19, 579-606, 2023 onset is too abrupt to be caused by insolation changes, and thus we can interpret that change with more confidence.The scatter in the late Holocene WDC 86 Kr xs data exceeds the stated analytical precision.Potential explanations include (1) an underestimation of the true analytical precision, (2) interannual to decadal variations in storminess at WDC, and (3) aliasing of centimeter-scale variations in ice core 86 Kr xs linked to layering in firn microstructural properties.Understanding the cause of this relatively high scatter in the 86 Kr xs records will require more work, in particular measurements of several high-resolution 86 Kr xs records in various sectors of Antarctica.

Barometric variability at WAIS Divide during the last deglaciation
In the present day, synoptic-scale pressure variability at WAIS Divide is correlated with zonal wind strength along the southern margin of the SPJ (Sect.4).In our interpretation, a more negative 86 Kr xs reflects a strengthening or southward shift of the SPJ in the Pacific sector.Here we provide a climatic interpretation of the deglacial WDC 86 Kr xs record and suggest that variations in synoptic variability at WDC are linked to movement of the ITCZ on millennial and orbital timescales.
The main features of the deglacial WDC 86 Kr xs record listed in Sect.5.1 resemble similar features seen in records of (sub)tropical hydrology and monsoon strength, such as the speleothem calcite δ 18 O records from Hulu Cave, China (Fig. 7c), and from Botuvera cave, southern Brazil (Fig. 7d), which are thought to reflect the intensity of the East Asian and South American summer monsoons, respectively (Cruz et al., 2005;Wang et al., 2001Wang et al., , 2007)).These two monsoon records are anticorrelated, showing opposing rainfall trends between the NH and SH on both orbital and millennial timescales.This pattern is commonly attributed to displacement of the mean meridional position of the ITCZ (Chiang and Friedman, 2012;McGee et al., 2014;Schneider et al., 2014), driven by hemispheric temperature differences (Fig. 7b).On orbital timescales such ITCZ migration has a strong precessional component, moving towards the hemisphere with more intense summer peak insolation; on millennial timescales the ITCZ responds to abrupt North Atlantic climate change associated with the D-O and Heinrich cycles (Broccoli et al., 2006;Chiang and Bitz, 2005;Wang et al., 2001), which are in turn linked to changes in meridional heat transport by the Atlantic meridional overturning circulation, or AMOC (Lynch-Stieglitz, 2017;Rahmstorf, 2002).
Changes in mean ITCZ position have a strong influence on the structure and strength of the SH jets.During periods when the NH is relatively cold (such as D-O stadials or periods with negative orbital precession index) the ITCZ is displaced southward and the SH Hadley cell is weakened, thereby also weakening the SH upper-tropospheric subtropical jet (Ceppi et al., 2013;Chiang et al., 2014).The re-verse is also true, with the ITCZ shifted northward during NH warmth, associated with a strengthening of the SH Hadley cell and STJ.In a range of model simulations (Ceppi et al., 2013;Lee and Kim, 2003;Lee et al., 2011;Pedro et al., 2018) the weakening of the SH STJ (as during NH cold) is furthermore accompanied by a strengthening and/or southward shift of the SPJ as well as the eddy-driven jet and SH westerly winds.Recently, ice core observations have confirmed that in-phase shifts in the position of the SHWs occur during the D-O cycle in parallel to those of the ITCZ (Buizert et al., 2018;Markle et al., 2017).Marine records of fluvial sediment runoff off the Chilean coast suggest precessionphased movement of the South Pacific SPJ, again in parallel to the ITCZ movement (Lamy et al., 2019).
The SAM index reflects the meridional position of the SHWs and eddy-driven jet.During positive SAM phases the SHWs are displaced poleward and during negative phases equatorward.Present-day month-to-month changes in the SAM index represent a mode of internal variability, with anomalies persisting for only weeks to months -the timescale is longest in late spring and early summer, reflecting a stronger planetary wave-mean flow interaction (Simpson et al., 2011;Thompson and Wallace, 2000).By contrast, shifts in the ITCZ and SH jet structure on millennial and orbital timescales have a much longer lifetime and different dynamics, being driven from the tropics via hemispherically asymmetric changes in Hadley cell and STJ strength.Therefore, present-day SAM internal variability is not expected to be a good analog for past changes in SHW position.We find that the present-day SAM month-to-month internal variability mainly impacts synoptic variability over the Southern Ocean and does not have a statistically significant impact at WDC (Table 2).Such variability is likely to have also occurred during other climatic regimes, possibly just centered around a mean SHW position that is displaced meridionally relative to today.At first glance it may appear contradictory to state, as we do, that synoptic activity at WDC is not sensitive to the SAM while also suggesting that during the last deglaciation synoptic activity at WDC is linked to changes in the position of the SH eddy-driven jet and westerlies.Based on the considerations above, both claims may be true without contradiction.
Besides secular changes to the SPJ position and strength linked to meridional ITCZ movement, WDC 86 Kr xs may also have imprints from ENSO and tropical Pacific climate.Our analysis suggests a weak but statistically significant link to common ENSO indicators (Table 2).Increased synoptic activity at WDC is linked to enhanced convection in the central and eastern tropical Pacific, which may be due to enhanced frequency or intensity of El Niño events or a mean climate state that is more El Niño-like; it seems likely that the Pacific mean state and ENSO variability are strongly linked (Salau et al., 2012), and the distinction may be irrelevant.
The key features of the WDC 86 Kr xs record are compatible with paleo-ENSO changes commonly described in the liter-ature.A majority of Holocene ENSO reconstructions (Conroy et al., 2008;Driscoll et al., 2014;Koutavas et al., 2006;Moy et al., 2002;Riedinger et al., 2002;Sadekov et al., 2013) and a wide range of climate model simulations (Braconnot et al., 2012;Cane, 2005;Clement et al., 2000;Liu et al., 2000Liu et al., , 2014;;Zheng et al., 2008) all suggest weakened ENSO activity during the early and mid-Holocene, a time with reduced WDC synoptic activity.For example, Fig. 7f shows the number of El Niño events per century (with trend line) reconstructed from inorganic clastic laminae in sediments from Laguna Pallcacocha, Ecuador, a region strongly affected by ENSO (Moy et al., 2002).Likewise, it has been suggested that the SST gradient between the western Pacific warm pool and eastern Pacific cold tongue was enhanced during the mid-Holocene, perhaps indicating a more La Niña-like mean climate state (Koutavas et al., 2002;Sadekov et al., 2013).
Going from the early Holocene to the Younger Dryas (YD), we observe a large increase in WDC synoptic activity.Enhanced ENSO activity during Heinrich stadials is generally supported by climate model simulations (Braconnot et al., 2012;Merkel et al., 2010;Timmermann et al., 2007) and by limited proxy evidence for stadial periods more broadly (Stott et al., 2002).Enhanced ENSO variability during the deglaciation is also found by Sadekov et al. (2013), although their record lacks the temporal resolution to resolve the individual stages.The zonal SST gradient in the equatorial Pacific further reaches a minimum during HS1, also consistent with higher El Niño intensity (Sadekov et al., 2013).
The observed variations in 86 Kr xs and implied changes in WDC synoptic activity may thus have two contributions: (1) ITCZ-driven changes to the South Pacific SPJ position and (2) changes to ENSO activity.Based on previous work, we argue these two amplify one another in driving WDC storminess, yet we expect the former to make the larger contribution.To disentangle zonally uniform changes to the SPJ from changes specific to the Pacific sector (such as ENSO and the split jet), 86 Kr xs records from different sectors of Antarctica are needed.Replication of the deglacial and Holocene WDC 86 Kr xs record presented here is also a high priority at both the WDC itself and the nearby SDM and RICE cores to validate that the signals we describe and interpret here are indeed real and regional in scale.
The position of the SHWs during the LGM has been a topic of much scientific inquiry.Proxy data have been interpreted to show a northward LGM shift of the SHWswith other scenarios, including no change at all, not excluded by the data (Kohfeld et al., 2013).Such a shift is not supported by most climate models (Rojas et al., 2009;Sime et al., 2013).Our 86 Kr xs record suggests LGM synoptic activity in West Antarctica to be comparable to today after accounting for site elevation effects (the elevation effect on 86 Kr xs is within the analytical error).This would be consistent with a Pacific SPJ position similar to today.Note that our site is mostly sensitive to the position of the southern edge of the SPJ and cannot meaningfully constrain changes to the seasonality, width, and/or northern edge of the storm tracks.Therefore, it is not a priori clear whether our observations can be extrapolated to more general statements about SHW position and strength during the LGM.Our data suggest that SPJ movement follows insolation and the ITCZ position, and hence the LGM period may not be a good target for studying SHW movement in the first place given that it has a precession index similar to the present day.
Changes to the SPJ and its associated westerly surface winds have implications for ocean circulation and marine productivity in the Southern Ocean via wind-driven upwelling.Opal flux records from the Antarctic zone (Fig. 7g), reflecting diatom productivity, are commonly interpreted as a proxy for such upwelling -with enhanced upwelling during southward displacement of the SHWs (Anderson et al., 2009).Here we only show records from the Pacific sector, given that we find WDC 86 Kr xs to reflect purely local SPJ dynamics (Fig. 4a).Both published records suggest enhanced upwelling during the deglaciation (Fig. 7g), consistent with a southward-shifted Pacific SPJ and enhanced storminess at WDC.The record from core PS75/072-4 (blue curve) further indicates an increasing productivity trend through the Holocene (Studer et al., 2018), which is accompanied by a rise in surface nitrogen availability (reconstructed from diatom-bound nitrogen isotopic composition, not shown); this Holocene trend matches our finding of increasing WDC storminess and, by inference, an increasingly southern position of the Pacific SPJ and SHWs.We thus conclude that our interpretation of WDC 86 Kr xs reflecting SPJ movement in parallel with the ITCZ is broadly consistent with indicators of wind-driven upwelling in the Pacific Antarctic zone.

Conclusions
Here we present a new gas-phase ice core climate proxy, Kr-86 excess, that reflects time-averaged surface pressure variability at the site driven by synoptic activity.Surface pressure variability weakly disturbs the gravitational settling and enrichment of the noble gas isotope ratios δ 86 Kr and δ 40 Ar via barometric pumping.Owing to its higher diffusion coefficient, argon is less affected by this process than krypton is, and therefore the difference δ 86 Kr-δ 40 Ar is a measure of synoptic activity.
This interpretation is supported by a calibration study in which we measure 86 Kr xs in late Holocene ice core samples from 11 Antarctic ice cores and 1 Greenland ice core that represent a wide range of synoptic activity in the modern climate.Two of the Antarctic cores were rejected due to clear evidence of refrozen meltwater.We find a strong correlation (r = −0.94when using site mean data and r = −0.83when using individual samples, p < 0.01) between ice core 86  Using atmospheric reanalysis data, we show that synopticscale barometric variability in Antarctica is primarily linked to the position and/or strength of the southern edge of the eddy-driven subpolar jet (SPJ, also called the polar front jet) with a southward SPJ displacement enhancing synopticscale surface pressure variability in Antarctica.The commonly defined modes of large-scale atmospheric variability, such as the southern annular mode and the Pacific-South American Pattern, impact Antarctic only weakly as they are weighted towards the midlatitudes; the exception is the Antarctic Peninsula, where synoptic activity is wellcorrelated with the southern annular mode (r = 0.68).Sites in the Amundsen and Ross Sea sectors are weakly linked to tropical Pacific climate and ENSO (r = 0.31 to r = 0.43).
We present a new record of 86 Kr xs from the WAIS Divide ice core in West Antarctica that covers the last 24 kyr including the LGM, deglaciation, and Holocene.West Antarctic synoptic activity is slightly below modern levels during the Last Glacial Maximum (LGM), increases during the Heinrich Stadial 1 and Younger Dryas North Atlantic cold period, weakens abruptly at the Holocene onset, remains low during the early and mid-Holocene (up to ∼ 17 % below modern), and gradually increases to its modern value.The WDC 86 Kr xs record resembles records of tropical hydrology and monsoon intensity that are commonly thought to reflect the meridional position of the ITCZ; the sense of the correlation is that WDC synoptic activity is weak when the ITCZ is in its northward position, and vice versa.We interpret the record to reflect migrations of the eddy-driven SPJ in parallel with those of the ITCZ (Ceppi et al., 2013).Secondary influences may come from tropical Pacific climate and ENSO activity.Our 86 Kr xs record is consistent with weakened ENSO activity (or a more La Niña-like mean state) during the middle and early Holocene and enhanced ENSO activity during NH stadial periods -both these features have been described in the paleo-ENSO literature.The inferred changes to the SPJ are broadly consistent with proxies that indicate enhanced wind-driven upwelling in the Pacific Antarctic zone during NH cold stadial periods.
Kr-86 excess is a new and potentially useful ice core proxy with the ability to enhance our understanding of past atmospheric circulation.More work to better understand this proxy is warranted, and presently the conclusions of this paper should be considered tentative.In particular, replication of the deglacial Kr-86 excess record presented here in nearby cores is needed before these results can be interpreted with confidence.A full list of suggested follow-up studies is given in Sect.3.3.Despite the many challenges of Kr-86 excess, its further development is worthwhile owing to the dearth of available proxies for reconstructing SH extratropical atmospheric circulation.

A1 Gas loss correction
Gas loss processes artificially enrich the δ 40 Ar isotopic ratio used to calculate 86 Kr xs (Kobashi et al., 2008b;Severinghaus et al., 2003Severinghaus et al., , 2009)).Figure A1b shows the relationships between the two most common gas loss proxies, δO 2 /N 2 and δAr/N 2 , for all samples in the calibration dataset.We find a slope close to the 2 : 1 slope commonly reported in the literature (Bender et al., 1995); the exception is the DE08-OH site where the data fall on a 1 : 1 slope.Depletion in fugitive gases (such as O 2 and Ar) represents the sum of losses during bubble closure in the firn (Bender, 2002;Huber et al., 2006;Severinghaus and Battle, 2006) and those during drilling, handling, storage, and analysis of the samples (Ikeda-Fukazawa et al., 2005).The patterns are inconsistent with storage conditions alone -for example, the DF and EDC cores were stored very cold and SP drilled very recently; yet all three have strong δO 2 /N 2 and δAr/N 2 depletion.Natural gas loss from the firn and artifactual loss during drilling likely dominate the signal.The DE08-OH samples were drydrilled and suffered from poor ice quality for the most depleted samples, which may explain the alternate 1 : 1 slope at the site (Appendix B); note though that a recent work suggests a ∼ 5 : 1 slope for post-coring gas loss (Oyabu et al., 2021).The DE08-OH samples were also analyzed differently from those at other sites, with δO 2 /N 2 and δAr/N 2 measurements performed on a separate smaller ice piece (see Sect. 2.2); the greater surface-to-volume ratio of such small samples may result in greater gas fractionation while evacuating the sample flasks in the laboratory.Severinghaus et al. (2009) hypothesize that the apparent 2 : 1 slope of δO 2 /N 2 to δAr/N 2 depletion is a combination of two mechanisms: size-dependent fractionation during diffusion through the ice lattice and mass-dependent fractionation (such as molecular or Knudsen diffusion) within ice fractures.In this interpretation, the exact slope would depend on the relative contribution of each process to the total gas loss.It is improbable that both processes would occur in the same ratio at such a wide variety of sites; the 2 : 1 slope is thus more likely an attribute of the gas diffusion rate of gases through ice itself, which is strongly size-dependent and weakly mass-dependent (Battle et al., 2011).
Gas loss is well-known to enrich ice samples in δ 18 O-O 2 , and following Severinghaus et al. (2009) we plot δ 18 O (corrected for gravity and small atmospheric δ 18 O atm variations) against gravitationally corrected δO 2 /N 2 in Fig. A1c.We find a slope of 3.5 per meg enrichment in δ 18 O per ‰ of δO 2 /N 2 gas loss.This is less than values reported elsewhere (Severinghaus et al., 2009) but provides further evidence for mass-dependent fractionation during gas loss.Our core-top dataset further suggests a correlation between δ 40 Ar − 4 × δ 15 N (a measure of δ 40 Ar enrichment impacted by both thermal fractionation and gas loss) and gravitationally corrected δAr/N 2 (Fig. A1d), suggesting Ar loss leads to enrichment of the remaining δ 40 Ar.
The value of ε 40 = −0.008agrees reasonably well with other studies.Kobashi et al. (2008b) compare replicate sample pairs to back-out gas loss and find (statistically significant) correlations between δ 40 Ar enrichment and δAr/N 2 (again, which is similar to δO 2 /N 2 -δAr/N 2 ).Kobashi et al. (2008b) find ε 40 values of −0.006, −0.005, and +0.007, depending on the depth range and analytical campaign evaluated.The positive value is surprising given that most observations, as well as theory, suggest ε 40 should be negative -we consider this a spurious result given the weak δ 40 Ar-δAr/N 2 correlation in that particular dataset.the Siple Dome ice core (Severinghaus et al., 2003), regressing δ 40 Ar against δKr/Ar gives a slope of +0.007; this implies ε 40 = −0.007, in good agreement with our findings.Last, our core-top data suggest δ 40 Ar enrichment with an ε 40 value of −0.0072 (Fig. A1d), also in good agreement with Byrd.Given the uncertainty in the gas loss parameter, we verify that our results are valid for a wide range of ε 40 values (Fig. 3b).
The downcore WDC 86 Kr xs was measured over five separate measurement campaigns, with δO 2 /N 2 and δAr/N 2 data available only for the first three campaigns.
In order to provide a consistent gas loss correction to the five measurement campaigns, including campaigns 4 and 5 for which no δO 2 /N 2 or δAr/N 2 data are available, we fit a third-order polynomial to all available gravitationally corrected WDC δO 2 /N 2 -δAr/N 2 data (Fig. A5a).We can then calculate the expected WDC δO 2 /N 2 -δAr/N 2 at any given age, also in the absence of δO 2 /N 2 and δAr/N 2 data.For consistency, we use this correction method for all data seen in Fig. 6.Note that the WDC δO 2 /N 2 -δAr/N 2 values are small for all ages and that the gas loss correction is therefore small for this site.

A2 Thermal correction
In the presence of a temperature gradient, thermal diffusion causes isotopic enrichment towards the colder location.The thermal diffusion sensitivity in units of ‰ K −1 for the various gases is given as (Grachev and Severinghaus, 2003a, b;Kawamura et al., 2013) follows.The samples from the calibration dataset are from the climatically stable late Holocene period and typically close together in depth; the uncertainty in the T estimation for individual samples therefore exceeds the temporal variability in T .
To reduce the uncertainty in the thermal correction we estimate T for individual samples using Eq.(A2) and for each site average the available data to get a site average firn temperature gradient T .The thermal correction is then given by the following.The two methods are compared in Figs.A3c (individual sample T ) and Fig. A3d (site mean T ); it is clear that the T approach reduces the spread in 86 Kr xs (error bars), but not its mean (white dots).The T estimates in individual samples are subject to errors in the isotopic measurements; some of these errors will cancel out in the T .
For the downcore WDC record through the deglaciation we can no longer assume a stationary T ; we instead rely on dynamic firn densification model simulations of T (Buizert et al., 2015).A comparison of the simulated and databased T is shown in Fig. A5 for WDC.The data clearly show a lot more scatter and variability than the simulations do.We interpret this mainly as analytical noise in the δ 15 N

A3 Elevation correction
To correct the deglacial WAIS Divide record for elevation changes, here we estimate the 86 Kr xs dependence on site elevation using the calibration dataset.Note that elevation and synoptic activity are strongly correlated for the investigated sites (r = −0.86),with synoptic activity decreasing with elevation because the cyclonic systems do not penetrate deeply into the Antarctic interior.Figure A6 shows the result of this exercise.We find a slope of 34 per meg ‰ −1 of 86 Kr xs per 1000 m of elevation change, with a correlation of r = 0.96 when considering site mean 86 Kr xs and r = 0.86 when considering individual samples.Note that the GISP2 site is not included in the analysis because it is in Greenland where the elevation-86 Kr xs relationship may be different from Antarctica -it does, however, fit the Antarctic trend rather well.We further use the simulated WAIS Divide elevation history (Golledge et al., 2014), which simulates an LGM elevation of around 300 m higher than at present at WAIS Divide.The Law Dome DE08-OH site is a revisit of the DE08 site, drilled in the 2018-2019 austral summer Antarctic field season.We have samples from two separate cores: (1) 13 24 cm long samples from a 10 cm diameter core going from 97 to 193 m depth at ∼ 8 m sample spacing and (2) 8 6 cm long samples from a 24 cm diameter core going from 97.6 to 99.8 m depth at 30 cm sample spacing.The purpose of the first set was to determine possible long-term variations in 86 Kr xs ; the purpose of the second set was to assess whether there are sub-annual variations in 86 Kr xs due to the seasonality in firn properties and bubble trapping.
Both cores were dry-drilled (i.e., no drill liquid was used).The 10 cm diameter core used was drilled at the beginning of the field season and the 24 cm diameter core at the end of the field season.Prior to shipment off the continent, both cores were stored in a chest freezer at Casey Station; due to a miscommunication this freezer was set to −20 • C rather than −26 • C, yet the ice is believed to have stayed below −18 • C.
Both DE08-OH cores experienced more gas loss than the original DE08 core that we also sampled (Fig. A1b).In particular, the samples from the 10 cm diameter core were strongly depleted in δAr/N 2 , with the most extreme gas loss seen for the deepest samples for which the ice quality was poorest.
Figure B1 shows the high-resolution sub-annual DE08-OH sampling.The data were corrected for gas loss and thermal fractionation using a site mean temperature gradient of T = −1.6 • C, possibly related to a rectifier effect (Morgan et al., 2022).We find strong (5-fold) variations in 86  sub-annual timescales.With an expected annual layer thickness of around 1.3 m at this depth, it appears as though there may be an annual-scale variation in 86 Kr xs ; the dataset has insufficient length to establish this firmly.We refrain from interpreting the long-term variations in 86 Kr xs in the 10 cm diameter core for two reasons.First, given the strong sub-annual variations seen in the highresolution sampling, it is unavoidable that we are aliasing the underlying signal in the core.Second, the 10 cm diameter core suffers from strong gas loss (depleted δAr/N 2 ).We attribute this primarily to the dry drilling and imperfect sample storage conditions.Perhaps the greater stresses during drilling a 10 cm core (compared to the 24 cm diameter core) result in more microfractures and gas loss.

Figure 1 .
Figure 1.Idealized firn air transport model experiments of 86 Kr xs .Firn density is calculated usingHerron and Langway (1980) and the diffusivity usingSchwander (1989).(a) Simulations using a fraction of dispersive mixing of f = 0 (left panel), f = 0.1 (middle panel), and f = 0.2 (right panel) for a hypothetical site with an accumulation rate of A = 2 cm a −1 ice equivalent and mean annual temperature T = −60 • C. At dispersive fraction f , effective molecular diffusivity of all gases is multiplied by (1-f ) and dispersive mixing for all gases is set equal to f times the effective molecular diffusivity of CO 2 .(b) Isotopic disequilibrium as a function of dispersive mixing intensity at two different firn thicknesses of around 100 m (dashed, A = 2 cm a −1 and T = −60 • C) and 50 m (solid, A = 2 cm a −1 and T = −43 • C).We compare isotopic disequilibrium without (blue, left axis) and with (orange, right axis) normalization.(c) Simulations at 10 % dispersive mixing, with each dot representing different climatic conditions.The accumulation rate is A = 2 cm a −1 ice equivalent, and mean annual temperature is changed from −60 to −30 • C in steps of 5 • C.

Figure 3 .
Figure 3. Calibrating Kr-86 excess.(a) 86 Kr xs as a function of for the calibration dataset.Circles give the site mean, and the error bars denote the ±1σ standard deviation between samples (uncertainty in corrections and measurements not included).The number of samples at each site is given in Table 1.The Pearson correlation coefficient is r = −0.94when considering site data means and r = −0.83when considering all individual samples.Data are corrected for gas loss using ε 40 = −0.008(Appendix A1) and corrected for thermal fractionation using site mean 15 N excess (Appendix A2).The calibration curve for 86 Kr xs15 is identical in this case, with slightly larger error bars.(b) Correlation of the calibration curve as a function of the gas loss correction scaling parameter ε 40 .The solid line gives the correlation for both site mean 86 Kr xs15 and 86 Kr xs40 (identical); the dashed lines show the correlation using individual samples for 86 K xs40 (blue) and 86 Kr xs15 (orange).Triangles denote the ε 40 estimate from the Byrd, Siple, and GISP2 ice cores (Fig.A2;Kobashi et al., 2008a;Severinghaus et al., 2003).

Figure 4 .
Figure 4. Zonal wind speed at 850 hPa (color shading, see scale bar) and 200 hPa (2 m s −1 contours) regressed onto surface synoptic activity at the Antarctic ice core sites of (a) WAIS Divide, (b) Law Dome (DE08, DE08-OH, and DSSW20K), (c) Dome Fuji, and (d) James Ross Island.Yellow dots mark the ice core locations.

Figure 5 .
Figure 5. Modes of SH extratropical atmospheric variability and their link to synoptic-scale surface pressure variability in Antarctica.(a) Annual mean in units of percent per day (% d −1 ).(b) Colors show correlation between and the southern annular mode (SAM) index, with the 500 hPa geopotential height anomalies superimposed in 10 m contours (positive contours solid, negative contours dashed).(c) As panel (b), but for the Pacific-South American Pattern 1 (PSA1).(d) As panel (b), but for the Pacific-South American Pattern 2 (PSA2).SAM, PSA1, and PSA2 are defined as the first, second, and third EOFs (empirical orthogonal functions), respectively, of the 500 hPa geopotential height anomalies in 20-90 • S monthly values in the 1979-2017 ERA-Interim reanalysis (Deeet al., 2011).In all panels the latitude of maximum is denoted by the green line.

Figure 6 .
Figure 6.WAIS Divide Kr-86 excess records through the last deglaciation.(a) WDC 86 Kr xs40 data from the five measurement campaigns.The gray curve shows a Gaussian smoothing curve to the combined data from the first two campaigns; the light gray shaded area shows the ±1σ uncertainty envelope based on a 10 000-iteration Monte Carlo sampling of the errors and uncertainties.The WDC calibration data are shown as gray circles for comparison.(b) As in panel (a), but for 86 Kr xs15 .For campaigns 4 and 5 the sample was not split, and no δ 15 N data are available.The Heinrich Stadial 1 and Younger Dryas North Atlantic cold periods are marked in yellow.Thermal corrections in the WDC 86 Kr xs records are based on firn model simulations; gas loss correction is based on a third-order polynomial fit to the WDC gravity-corrected δO 2 /N 2 − −δAr/N 2 (Fig. A5).

Figure 7 .
Figure 7. Climate records through the last deglaciation with the Heinrich Stadial 1 (HS1) and Younger Dryas (YD) North Atlantic cold periods marked in yellow.(a) Greenland Summit ice core stable water isotope ratio δ 18 O, here the average of the GISP2 and GRIP ice cores(Grootes et al., 1993).(b) Hemispheric temperature difference(McGee et al., 2014) based on global proxy compilations for the Holocene(Marcott et al., 2013) and last deglaciation(Shakun et al., 2012).(c) Speleothem calcite δ 18 O from Hulu and Dongge caves, China, as a proxy for East Asian summer monsoon strength(Dykoski et al., 2005;Wang et al., 2001).Superimposed is summer solstice (21 June) insolation at 30 • N. (d) Speleothem calcite δ 18 O from Botuvera cave, southern Brazil, as a proxy for South American summer monsoon strength(Cruz et al., 2005;Wang et al., 2007).(e) Kr-86 excess record from WAIS Divide (this study), corrected for gas loss and thermal fractionation (Appendix A).The center line and shaded envelope show the mean and ±1σ uncertainty interval of a 10 000-iteration Monte Carlo smoothing exercise (see text).The dotted red line equals the center line with a correction for elevation change applied (Appendix A) using a simulated elevation history(Golledge et al., 2014).(f) Number of El Niño events per century from laminations in sediments from Laguna Pallcacocha, Ecuador(Moy et al., 2002).(g) Th-normalized opal flux in the Pacific Antarctic zone (south of the polar front) from cores NBP9802-6PC1 (turquoise; 169.98 • W, 61.88 • S) and PS75/072-4 (blue; 151.22 • W, 57.56 • S), reflecting local productivity and (wind-driven) upwelling(Chase et al., 2003;Studer et al., 2015).All water isotope data in this figure are on the Vienna Standard Mean Ocean Water (V-SMOW) scale.Arrows show the direction of increased monsoon strength and synoptic activity.

Figure A1 .
Figure A1.Elemental ratios in the 11-site calibration study of late Holocene samples.(a) δXe/N 2 vs. δ 40 Ar in all ice core samples.δ 40 Ar is used solely to illustrate gravitational enrichment, and a similar picture arises when plotted against any isotopic pair.Refrozen meltwater (elevated δXe/N 2 ) was seen in all samples from the Antarctic Peninsula (James Ross Island and Bruce Plateau sites), despite selecting samples free of visible melt features.(b) The relationship between the commonly used gas loss proxies δO 2 /N 2 and δAr/N 2 corrected for gravity.(c) Enrichment in δ 18 O (corrected for gravity and atmospheric δ 18 O atm ) plotted against gravity-corrected δO 2 /N 2 .(d) δ 40 Ar enrichment plotted against gravity-corrected δAr/N 2 .In all panels gravitational correction is applied by subtracting δ 15 N times the atomic mass unit difference.

Figure A2 .
Figure A2.Argon isotopic enrichment due to gas loss in the Byrd core used to determine the δ 40 Ar gas loss correction (Appendix A1).The enrichment in δ 40 Ar plotted as a function of gravitationally corrected (δO 2 /N 2 -δAr/N 2 ) measured in the deep Antarctic Byrd ice core, which suffered heavy gas loss.Ice samples were analyzed in the Bender Lab at the University of Rhode Island by Jeff Severinghaus in 1997.The slope of the least-square fit is ε 40 = −0.008.The data point in parentheses is treated as an outlier and excluded from the fitting.
thermal gradient T in the firn using N-15 excess(Severinghaus et al., 1998): the δ 40 Ar gas loss correction from Eq. (A1).Positive values of T indicate that the surface is warmer than the firn-ice transition.The T then in turn allows us to estimate the thermal corrections.

Figure A3 .
Figure A3.Influence of gas loss and thermal correction on the 86 Kr xs40 calibration.We plot 86 Kr xs40 as a function of (a) without any data corrections applied, (b) with only the gas loss correction applied (ε 40 = −0.008),(c) with only the thermal correction applied using individual sample T , (d) with only the thermal correction applied using individual site mean T , (e) with both gas loss and thermal corrections applied using individual sample T , and (f) with both gas loss and thermal corrections applied using site mean T .In each panel the correlations with are listed for the site average and individual sample with the latter in parentheses.For all correlations p < 0.05.

Figure A4 .
Figure A4.Same as Fig. A3, but for 86 Kr xs15 .Note that the gas loss correction (b) does not impact 86 Kr xs15 .For all correlations p < 0.05, except for panels (a) and (b) for which p = 0.16 for the site average correlation.

Figure A5 .
Figure A5.Gas loss and thermal corrections for the WDC time series.(a) WDC gravity-corrected δO 2 /N 2 -δAr/N 2 as a measure of gas loss.The solid line is a third-order polynomial fit to the data; the dashed lines give a ±1 ‰ range around the fit, which captures the majority of the data.(b) The T correction applied to the downcore records.The blue envelope shows the ±2σ range of thermal correction scenarios in the Monte Carlo sampling, together with the mean (blue line).Gray dots show WDC T estimates from available 15 N excess data, with the red curve being a Gaussian smoothing function to the data.

Figure A6 .
Figure A6.Kr-86 excess dependence on site elevation.The vertical axis is the 86 Kr xs .The linear fit has a slope of 34 per meg ‰ −1 per 1000 m elevation.

Figure B1 .
Figure B1.High-resolution sub-annual sampling of 86 Kr xs40 at the DE08-OH site.The annual layer thickness at this depth is around 1.3 m.

Table 1 .
Ice core sites used in this study, with N the number of samples included in the calibration study.See the main text for acronyms.
Kr xs and barometric variability at the site.Current limitations of the new 86 Kr xs proxy are the following: (1) it requires relatively large and nontrivial correchttps://doi.org/10.5194/cp-19-579-2023Clim.Past, 19, 579-606, 2023 tions for gas loss and thermal fractionation; (2) it is moderately sensitive to changes in convective zone thickness; (3) firn air transport models cannot simulate the magnitude of 86 Kr xs anomalies measured in ice samples; (4) firn air samples show smaller 86 Kr xs anomalies than ice samples from the same site do; (5) it may be sensitive to the degree of density layering at the site, as a comparison of the nearby Law The other two values of ε 40 are in reasonable agreement with the Byrd value.For https://doi.org/10.5194/cp-19-579-2023Clim.Past, 19, 579-606, 2023