1. Dry and wet extraction COS measurements
COS measurements were conducted at the UCI ice core trace gas laboratory using established methods60 to analyze 10-20 cm long samples from the designated gas cross-section of the SPC14 ice core. Air entrapped in the ice core samples is extracted using dry and wet extraction methods. In dry extraction, ice core air is liberated via mechanical shredding of frozen ice samples in four identical stainless-steel vacuum chambers (6 L), each fitted with a flat, sharp cutting surface. Samples are placed on the cutters inside the vacuum chambers and shredded by linear oscillations (15 cm throw at 1 Hz for 20 min) inside a -50°C freezer. In wet extraction, ice core air is liberated by melting the samples in two identical glass vessels at 30°C over 30-40 min. The SPC14 dry extraction samples were 400-500 g and wet extraction samples were 300-350 g. Melting releases close to 100% of the ice core air. The dry extraction efficiency is about 60% for bubbly ice samples and 40% for full clathrate ice, with a drop in extracted air amount across the bubble-clathrate transition zone (BCTZ) (EDFig. 1a).
Following either wet or dry extraction, the sample air is frozen in a stainless-steel tube dipped in liquid helium (4K) for transfer to the analysis inlet system. In the inlet system, trace gases in the ice core air are pre-concentrated on a glass-bead trap using liquid nitrogen (77K), then transferred onto a capillary trap (77K) before being thermally injected onto a DB-1 gas chromatography (GC) column. The GC is connected to a dual focusing magnetic sector mass spectrometer (MS) that operates at a mass resolution (m/Dm) of 8,000-10,000 at 5% peak height. The GC-MS system is calibrated with an isotope dilution technique using a 13C-labeled COS isotope standard (13COS) with ppt-level unlabeled 12COS standards. Both labeled and unlabeled ppt-level standards are prepared in our laboratory from ppb-level primary (long term) standards. The primary standards are also prepared in our laboratory from commercially available pure gases. A known amount of ppt-level 13COS internal standard is mixed with every sample during the pre-concentration step. Dry extraction-based COS measurements from other ice cores have been published19-20. The wet extraction system and methods have been developed for non-methane hydrocarbon measurements61. The SPC14 measurements presented here are the first COS data from the wet extraction system. All dry and wet extraction measurements are corrected for 2% gravitational enrichment based on hydrostatic equilibrium at the present-day firn bottom depth of 120 m at the South Pole.
There are no dry measurements from the initially anticipated BCTZ zone (900-1200 m) (Fig. 1a) to avoid measurement artifacts due to the low extraction efficiency of dry extraction20. The wet extraction method is more than 99% efficient in extracting ice core air, thus is not prone to artifacts in the BCTZ. However, COS is more soluble in water than air, resulting in a negative bias in wet measurements61 (Fig. 1a). We calculate the solubility correction empirically for each sample and correct the wet measurements (Methods-2). The solubility-corrected wet measurements exceed the dry measurements in the 800-900 m range, indicating a depletion in the dry measurements up to 100 m above the initially anticipated depth range of the BCTZ. We interpret this apparent depletion in dry measurements to be indicative of the onset of the BCTZ at about 800 m.
2. The solubility correction for the wet extraction measurements
The ratio of total moles of COS to what is left in the ice core melt water during wet extraction can be calculated by61:
In equation-1, mt is the total moles of COS, ml is moles of COS left dissolved in the melt water, Vg is the volume of the gas (headspace) in the melt chamber, Vl is the volume of the melt water, and α is the dimensionless effective solubility. We assume α to be the same for every sample melt. Vg can be substituted by volume of the melt chamber (Vc) minus Vl, and equation-1 can be rearranged into expressing mt/ml as a linear function of a single variable Vl and two constants (C1 and C2), given that Vc is also constant:
$$\frac{{m}_{t}}{{m}_{l}}={C}_{1}+{C}_{2}\frac{1}{{V}_{l}}$$
2
On the left side of equation-2, mt scales with the dry extraction COS mixing ratio (COSdry) and ml scales with the difference between the dry and wet (COSwet) extraction mixing ratios at the same depth, assuming solubility of air in the melt water is negligible compared to that of COS. Nicewonger et al.61 estimated 0.7% of the air could be left in the melt water for an average size sample at saturation equilibrium. The dry extraction mixing ratio at a given depth is estimated by linear interpolation from the smoothed dry extraction record (EDFig. 2). The Vl on the right hand side of equation-2 scales with the weight of the ice samples (Wts), assuming sample-to-sample density variations are negligible. Equation-2 can be rearranged to:
$$\frac{{COS}_{dry}}{{(COS}_{dry}-{COS}_{wet})}={C}_{3}+\frac{{C}_{4}}{{Wt}_{s}}$$
3
C 3 and C4 are constants that can be estimated from the linear regression of the inverse of Wts versus the left hand side of equation-3 using data from above 720 m (EDFig. 3a). We use a Bayesian errors-in-variables linear regression to estimate the solubility correction for the SPC14 wet extraction data due to nonuniform errors, although a least squares linear regression yields similar results (EDFig. 3). The errors-in-variables regression is conducted within a hierarchical Bayesian framework using the Stan probabilistic software package (mc-stan.org), in which all uncertain parameters are interpreted probabilistically62. Stan uses efficient Hamiltonian Markov Chain Monte Carlo (HMC) sampling. The Bayesian model assumes distributions for the x variable (1/Wts) using equation-4,
which describes the measured x as normally distributed with a mean equal to the true x and standard deviation xerr. The probability distribution of the linear relationship between x and y (left hand side of equation-3) variables is described by equation-5:
where parameters b and a are the unknown intercept (C3) and slope (C4) of the linear relationship. We do not impose any a priori limits on the parameters a and b. The true x parameter is required to be positive because sample weight cannot be less than zero. The posterior probability distributions are calculated via the Stan code (Supplementary Code) as executed in cmdStan 2.29 (https://github.com/stan-dev/cmdstan/releases). The HMC code uses 4 chains, 2000 iterations each for a total of 8000 evaluations of the posterior probability. The HMC diagnostic split \(\widehat{R}\) achieved \(\widehat{R}<1.01\) for the four chains, indicating convergence.
Once C3 and C4 are known, the solubility correction for wet extraction measurements (Fig. 1) can be estimated by inserting the COSwet and Wts data into equation-3 to calculate the solubility corrected COS mixing ratio represented by COSdry in equation-3. The calculations are carried out within the generated quantities block of the Stan code (Supplementary Code) using the full probability distributions of C3 and C4, and the uncertainty distributions of the x and y data. This allows for the propagation of all measurement and parametric uncertainty to the corrected wet extraction COS record. The average solubility correction (COSdry/COSwet) is a factor of 1.20 for samples above 720 m and 1.23 for samples below 720 m (EDFig. 3b).
3. Calculation of sea salt Na+ and non-sea salt Ca2+
We use established methods to infer the source contributions of Na+ and Ca2+ ions for the SPC14 ice core24. We partition Na+ and Ca2+ into sources from marine and crustal origin based on the respective Ca2+/Na+ mass ratios of 0.038 and 1.7863. Given these ratios and assuming that the sum of marine and terrestrial Na+ and Ca2+ equals what is measured in the laboratory, we use:
$${ssNa}^{+}=\frac{{Na}^{+}-\frac{{Ca}^{2+}}{{R}_{t}}}{1+\frac{{R}_{m}}{{R}_{t}}}$$
6
\({nssCa}^{2+}=\frac{{Ca}^{2+}-{R}_{m}{Na}^{+}}{1-\frac{{R}_{m}}{{R}_{t}}}\) (7)
where Rm = 0.038, Rt = 1.78, nssCa2+ is terrestrial non-sea salt Ca2+ and ssNa+ is sea salt Na+. ssNa and nssCa are commonly used as reliable proxies for ice core impurities of terrestrial and marine origin64.
4. ssNa as a proxy for COS production and the multiple analysis scenarios
Ice core impurities can lead to production of trace gases in the ice sheet22–23. When the accumulated production is comparable to the atmospheric levels of the measured gas, a correlation between the gas and the ice impurities will start to emerge. If the trace gas production from ice impurities occurs after bubble close-off, the excess gas should correlate with the impurity concentrations measured over the exact same depth range as the gas sample. Impurity induced gas production can also happen in the laboratory during the extraction of the ice core air, also resulting in correlation with the impurities measured over the exact same depth range as the gas sample. Trace gas production can also happen in the firn, in which case the excess trace gas is likely to be produced from impurities that cover a larger depth range than the exact length of the ice core samples because of vertical mixing of gases in the firn.
The dry COS measurements used in the production analysis are from the 1200–1751 m depth range (23–52 ky). These samples were about 20 cm long. The mean annual layer thickness in this depth range is 2 cm so each sample averages 10 years of ice chemistry. Linear regressions were conducted between COS and 10-y, 100-y, and 300-y averaged impurities (ssNa and nssCa) to determine whether marine or terrestrial impurities correlate with measured COS concentrations and to distinguish between potential production in the firn versus in the ice or the laboratory (EDFig. 4). In these analyses, we primarily rely on Bayesian errors-in-variables linear regressions (Methods-7) because both the x and y variables have non-uniform errors, although standard linear regressions were also conducted. The time averaged impurities are centered around the mid depths of the COS samples regardless of the averaging window.
We find significant slopes between COS and both ssNa and nssCa from the same depth regardless of the averaging window for the ice impurities, indicating some level of COS production from ice impurities. The slope is stronger and more significant with ssNa than nssCa for all averaging windows; we interpret this to be indicative of COS production being linked to marine sourced impurities. ssNa and nssCa are correlated in the SPC14 glacial ice, although the relationship is non-linear, resulting in the weaker COS-nssCa correlation. 10-y average ssNa yields a significant slope with COS, but the lowest R2 (EDFig. 4). The significance of the slope and R2 increase when the ssNa averaging window is increased to 100-y and slightly decline at 300-y. The higher frequency variability in ssNa does not correlate with the COS spikes in glacial ice, arguing against impurity driven COS production after bubble closure or during the extraction in the laboratory. We conclude that most of the COS production happens in the firn. The errors-in-variables slope of the 100-y ssNa vs. COS regression is used to correct the COS record and the results are presented as analysis scenario G1 (Figs. 2 and 3).
In three additional analysis scenarios (G2 through G4), we explore the uncertainty and robustness of the same-depth ssNa vs. COS relationship with respect to i) the COS spikes, ii) the COS measurement errors, and iii) the possible influence from a climate driven relationship between ssNa-COS. These additional analysis scenarios quantify the uncertainties arising from various data treatment decisions. In G2 (EDFig. 5d-f), we smooth the COS measurements with a 7-point moving average, which approximates to an 1143-y averaging window, and use the 7-point moving standard error of the measurements themselves as the y-variable error, in place of the original measurement errors. The ssNa measurements are also smoothed to 1100 years to match COS. Collectively, these choices eliminate the influence of the COS spikes from the regression via smoothing and introduce higher uncertainty to the y-variable where there is a higher occurrence of COS spikes. The R2 for G2 is significantly higher while the slope estimate increases slightly, indicating the same-depth ssNa-COS relationship is driven primarily by millennial scale variability and is not influenced by the higher frequency COS spikes.
The analysis scenario G3 (EDFig. 5g-i) differs from G2 primarily in that we reduce the weight of COS spikes on the regression by smoothing the COS record by a 7-point moving median. The treatment of y-variable errors also changes. We use a 7-point root mean square of the original COS measurement errors and introduce a constant multiplier (alpha) to allow the Bayesian inference to scale up the y-errors as needed (Methods-7). The G3 scenario demonstrates the impact of minimizing the influence of COS measurement uncertainties on the errors-in-variables analysis as the y-variable errors in G3 are more uniform by design and are scaled up by more than a factor of three, resulting in a closer slope estimate to the least squares linear regression of the data used in the G3 scenario.
In G4 (EDFig. 5j-l, 6d-i), we introduce a simultaneous regression of COS with ssNa from the same age (time). A same-age correlation between COS and ssNa would be indicative of a climate driven relationship. The delta-age is 1500–2700 years over the 1200–1751 m depth range (EDFig. 1b) corresponding to 30–54 m of ice (average annual layer thickness is 2 cm), suggesting only multi-millennial scale climate trends can hinder accurate quantification of the same-depth slope. We also conduct an independent same-age only regression of measured COS vs. 100-y smooth ssNa (EDFig. 6a-c, G1 for the same-age regression). The same-age regression slopes for G1 and G4 have mostly overlapping PDFs and similar means (EDFig. 6 c and f), revealing a negative correlation between ssNa and COS over multi-millennial time scales which impacts the slope of the same-depth relationship (EDFig. 5i). Detrending COS with the same-age ssNa-COS relationship increases the R2 of the same-depth regression to over 30% (Fig. 2a, EDFig. 5j). We consider the G4 results to be a more accurate estimate of the same-depth slope than G3 as it implicitly corrects for the influence of the climate driven variance in the record.
We conducted two similar data analysis scenarios (H1 and H2) for the dry extraction measurements from above 800 m. In H1 (EDFig. 7a-c, g-i), the same-depth and same-age correlations are identified independently using COS measurements (49-y resolution on average) and 100-y smooth ssNa. In H2 (EDFig. 7d-f, j-l), we apply a 7-point smoothing to the COS record (350 y smoothing window on average) and conduct simultaneous same-depth and same-age regressions vs. 300-y smooth ssNa. We find significant positive correlations between ssNa and COS from the same depth and from the same age. Quantification of the same-age relationship with the simultaneous regression does not significantly change the mean slope of the same-depth relationship but results in a narrower probability distribution.
We conducted two sensitivity tests for G1 and G4; first, to test that ssNa interpolated exactly to the COS measurements depths is indeed the best predictor of the COS production in the firn; second, to test the sensitivity of the regression results to random elimination of a subset of the measurements from the regression analysis. In the first test, we repeated the G1 analysis with 100-y smooth ssNa interpolated to 5 m and 10 m above and below the actual COS measurement depths (EDFig. 8a-f). The significance of the ssNa-COS relationship deteriorates with offsets in either direction, confirming that ssNa interpolated to the same depths as the COS measurements is the best predictor of the COS production, although this evidence alone is not sufficient to attribute the inferred production to a specific depth range in the firn.
In the second test, we repeated the G1 and G4 analyses two more times each with 20 subsets of the original COS data set (total of 40 simulations). Data in each subset is randomly picked from the original data without replacement to sample 80% and 60% of the data in the original data set. In G1 (EDFig. 8g), the mean slope value of the same-depth regression is not sensitive to the reduction in the amount of data included in the analysis but the shape of the distribution and the significance are. In G4 (EDFig. 8h,i), the simultaneous same-depth and same-age regression results are robust with respect to the reduction in the data amount, with some reduction in significance at 20% data reduction and additional reduction in significance at 40% data reduction, but with minimal impact on the mean value and the shape of the distribution.
5. Robustness of the same-age ssNa vs. COS regressions
We investigate whether correcting the COS record with the slope of the same-depth ssNa relationship results in a same-age correlation with ssNa in three analysis scenarios. In theory, this could explain the observed same-age negative correlation, in which case the same-age correlation cannot be interpreted as a climate signal. In the first scenario (EDFig. 9a), we interpolate an inverted (multiplied by -1 to simulate the ssNa correction to the COS record) 100-y smooth ssNa record to the SPC14 gas chronology at annual resolution, then regress it against the original 100-y ssNa record on ice chronology for the same depths over periods older than 23 ky, 25 ky, and 30 ky. Note that the dry extraction measurements from the glacial period extend from 23 ky to the bottom of the ice core. The second scenario (EDFig. 9b) is identical to the first except we used the 1100-y smooth ssNa record. The third scenario (EDFig. 9c) is a repeat of the second at lower and irregular resolution, using the 1100-y ssNa data corresponding to the COS measurement depths. All three tests reveal significant but small negative slopes for the 23 ky cut-off, but truncating the analyses with the 25 ky or 30 ky cut-offs results in either no significant correlation or a positive correlation.
These three tests can be compared with the same-age regressions between ssNa and COS. When the corrected COS from G4 is regressed with the 1100-y ssNa from the COS measurement depths, the slope is higher and it is more robust against the temporal cut-offs (EDFig. 9d). Repeats of the G4 regression itself using different temporal cut-offs reaffirm that the simultaneous same-age (EDFig. 9e) and same-depth (EDFig. 9f) regression results are robust with respect to the applied temporal cut-offs. Collectively, these tests show the same-age relationship between ssNa and COS prior to 23 ky result primarily from temporal trends that are inherent to the COS data set and do not result from the same-depth ssNa correction.
We note that the actual temporal trends in atmospheric COS might not closely track the mean of the posteriors of the corrections due to possible temporal variations in the slope, but the reported uncertainty ranges (Methods-7) account for a range of slopes within the posterior PDFs. The possible effects of more extreme temporal variability in the correction slopes can be inferred from the difference between the results for different analysis scenarios, with the caveat that part of the differences in slopes for the different analysis scenarios result from the reduction in dynamic ranges of the ssNa and COS records due to applied smoothing.
6. The COS rise concurrent with the deglaciation
We investigate whether there is any explicit statistical evidence that the increase in the inferred atmospheric COS record concurrent with the deglaciation during 19 − 10 ky is climate driven. As noted previously, the same-depth ssNa correction alone would result in a temporal COS trend that looks like the inverted ssNa record on the gas age scale (correction-ssNa) whereas the true climate signal in the ssNa record is evident on the ice chronology (climate-ssNa). The corrected COS record from G4 is regressed against the 1100-y smooth (G4) correction-ssNa and against the climate-ssNa (EDFig. 10). The R2 statistic for the regression versus the correction-ssNa is very high (EDFig. 10c), indicating the magnitude of the deglacial COS rise is set primarily by the ssNa correction. The R2 of the regression versus the climate-ssNa should be comparatively lower if the COS measurements merely introduce random noise into the corrected COS record, but the climate-R2 is higher than the correction-R2 (EDFig. 10c vs. d).
It is possible to estimate the probability that we would observe these climate-R2 and correction-R2 values coincidentally. To do so, we introduce random noise to the correction-ssNa, then regress against the actual correction-ssNa (itself without the noise), and separately against the climate-ssNa (also without noise) in 10,000 simulations, calculating the correction-R2 and climate-R2 for each instance. We find that the conditional probability of climate-R2 being higher than the correction-R2, given that climate-R2 is equal to or higher than 0.88, peaks at about p = 0.002 with 6 ppb (±1stdev) of added noise (EDFig. 10e). The probabilities are lower at lower noise because the correction-R2 approaches 1, and also lower at higher noise because the climate-R2 decreases rapidly. It is highly probable that the wet COS measurements from the deglaciation include a climate signal, which is retained during the same-depth ssNa correction, resulting in an inferred atmospheric COS record that is more similar to the deglacial climate signal than the inverse of the applied ssNa correction. This test does not rule out the possibility of temporal changes in the correction slope during the deglaciation, which is probable due to the large changes in the firn characteristics during this period21, but suggests that the effect of any such change does not negate the climate signal.
7. Bayesian methods for ssNa vs. COS regressions
All errors-in-variables linear regression analyses are conducted using a similar hierarchical Bayesian framework as the solubility correction for wet extraction data. The Bayesian model assumes distributions for the x variable (ssNa) based on equation-4. The probability distribution incorporating the linear relationship between x and y (COS) variables is characterized by equation-8, which is the same as equation-5 except it includes a third parameter α to scale the y-variable errors as needed depending on the analysis scenario.
Equations 4 and 8 are used for the regressions of x and y variable pairs for different cases during both the glacial period and the Holocene. We do not impose any a priori limits on the a, b, and α parameters. The true x (xtrue) parameter is required to be positive in the all regressions because ssNa cannot be less than zero. In three of the analysis scenarios (G1, G2, H1), the α parameter is constrained to a constant value 1, meaning the estimated y-variable uncertainties were not scaled in the Bayesian inference. The G3 regression includes the unconstrained α. The G4 and H2 regressions include simultaneous regressions of the same-depth and the same-age relationships as well as the unconstrained α. The simultaneous regression algorithm includes two instances of the equation-4 and two interdependent equations replace equation-5, characterizing the linear relationship for the same-depth and the same-age analyses:
The simultaneous regression includes 6 unknowns: x1 and x2 represent the true x from equation-4 for the same-depth and same-age regressions; a1, a2, b1, and b2 are slopes and intercepts for the same-depth and same-age regressions. The model accounts for the relationships in both spaces by estimating the parameters simultaneously. The simultaneous regression algorithm is the most complex Stan code used in our analyses (Supplementary Code). Simplified versions were used as needed for the different analysis scenarios described above. The regressions were conducted with cmdStan 2.29 (https://github.com/stan-dev/cmdstan/releases), using 4 chains, 2000 iterations each for a total of 8000 evaluations of the posterior probability. The HMC diagnostic split \(\widehat{R}\)(a measure of convergence between and within chains) of the slope posteriors for the G1, G2, G3, G4, H1, H2 scenarios achieved \(\widehat{R}\le 1.001\) with respective effective sample sizes (Neff ) of 4.9e3, 2.5e3, 7.1e3, 3.9e3, 3.9e3, 4.7e3. The convergence criteria are met65 when \(\widehat{R} <\) 1.01 at Neff > 400. The same diagnostics for all other Stan-based regressions also indicate convergence.
All slope PDFs shown in the manuscript display the posteriors of the slope parameter from the Bayesian inference. We also estimate a probability range for the corrected COS records for all regressions using the “generated quantities” block of the Stan code (Supplementary Code), accounting for the full range of uncertainties that arise from the slope PDFs and the x-variable ssNa. The 2σ uncertainty ranges shown in Fig. 3d for all scenarios represent ±2 stdev of the posterior probability distributions of the corrected COS records.
8. COS and ssNa errors for the Bayesian analyses
There are two primary factors that contribute to the error estimates of the individual COS measurements60. The first one is the uncertainty that arises from the reproducibility of calibration curves. A typical COS calibration curve for the GC-MS system is more uncertain at the low end of its dynamic range. This results in larger relative errors for lower COS mixing ratios and for smaller ice samples that yield less air for the analysis (EDFig. 1a). The second factor is the background COS that is inherent to the extraction vessels and vacuum lines used for the ice core air extractions. The background COS is characterized by clean N2 analyses through the ice core extraction system, which is subsequently subtracted from the ice core sample signal60. The variance in the background COS level introduces a larger relative uncertainty to smaller samples and lower mixing ratios much like the calibration curve uncertainty. As a result, the reported COS measurements have non-uniform errors (Fig. 3a).
ssNa was measured continuously in the SPC14 ice core, leading to much higher data density than the COS measurements24. We use a moving standard error of measurement means that fall within a 100-y window (higher or lower for other scenarios) to calculate the uncertainty in the ssNa concentration at the corresponding depth range in the ice core. The standard errors are calculated from unbiased (divided by N-1) standard deviations. There are on average 82 ssNa measurements per 100-y averaging window below 1200 m, which is the depth range corresponding to the glacial dry extraction measurements.
9. Comparison of the SPC14 COS record with COS records from other Antarctic ice cores
Comparison of the SPC14 record to previous ice core measurements is challenging. There is no other COS data set that does not require a hydrolysis loss correction or at a comparable resolution to the SPC14 record from the glacial period that would allow an investigation of impurity driven COS production impacting the measurements. We presume that some COS production occurs at every site and compare the SPC14 results with two other sets of ice core COS measurements that extend back into the last glacial period using the G4 correction slope estimate from the SPC14 ice core. These comparisons should be considered preliminary given that the general applicability of the same-depth ssNa-COS regression slopes has not been tested with data from the other ice core sites. A new high resolution data set, preferably from another cold East Antarctic site, is needed for a proper validation of the current results.
Taylor Dome M3C1 ice core
Measurements from the Taylor Dome M3C1 ice core were previously used to determine the kinetic parameters of in situ loss of COS19–20. During the Holocene, the loss corrected Taylor Dome data compares well with measured COS levels in the SPC14 ice core indicating the loss correction applied to the Taylor Dome record is reasonably accurate. However, the agreement is not good during the glacial period (EDFig. 11a). We apply a ssNa correction (using the SPC14 parameters) to the Taylor Dome measurements to explore the first order effects of a ssNa correction on the COS measurements (EDFig. 11b). The COS production occurs in the firn prior to significant loss, which is a slow process occurring over multi-millennial time scales19–20. Therefore the ssNa correction (G4) is applied to the loss-corrected COS record. Taylor Dome ssNa was calculated from the soluble Na and Ca measurements66 using equations-6 and 7.
After the ssNa correction, the Taylor Dome record agrees better with the SPC14 record during the glacial period from 45–20 ky and also suggests a steep increase during the deglaciation similar to the SPC14 record (EDFig. 11b). The timing of the deglacial COS increase is somewhat different for the Taylor Dome record versus the SPC14 record; the relatively minor discrepancies could result from multiple factors. It could be due to use of SPC14 specific ssNa correction parameters. The full range of uncertainty in the Taylor Dome loss correction arising from the uncertainties in the temperature histories have not been quantified19–20. The Taylor Dome chronology during the LGM and the deglaciation is susceptible to larger uncertainties than the SPC14 ice core, including known errors during the early Holocene because the current Taylor Dome gas chronology is tied to the outdated EDC-1 chronology67. Chronology errors impact the hydrolysis loss correction too as this correction is both temperature and time dependent.
WAIS Divide WDC-06A ice core
It was previously documented that the loss corrected WDC-06A record shows poor agreement with the Taylor Dome record during most of the Holocene except the last 1000 years (EDFig. 12a). This discrepancy was attributed primarily to the onset of BCTZ in the WDC-06A ice core; the Taylor Dome ice core is not deep enough to allow clathrate formation18. There is also a discrepancy between the WDC-06A and SPC14 measurements from 20–1 ky. Prior to 20 ky, the hydrolysis loss corrected WDC-06A record lies at the low end or just outside of the uncertainty band of the SPC14 measurements (EDFig. 12a). We apply the ssNa correction to the hydrolysis loss corrected WDC-06A record from below the presumed BCTZ using the SPC14 parameters (EDFig. 12b). The ssNa record for the WDC-06A is calculated using the total Na and Ca measurements68 in equations-6 and 7.
The ssNa corrected WDC-06A record agrees better with the SPC14 record during the glacial period and through 14 ky, capturing the first few thousand years of the COS rise during the glacial/interglacial transition. However, the WDC-06A does not show any indication of a peak from 20–18 ky. The WDC-06A record also displays a deep trough around 13 ky (EDFig. 12b). The SPC14 gas chronology is tied to the WDC-06A gas chronology by methane measurements and the relative age uncertainty between the two ice cores is 100 y or less for most of the core21. At least part of the discrepancy between the three ice cores during the late deglaciation is likely due to the fact that the SPC14 ssNa correction parameters are not broadly applicable to the Taylor Dome and WAIS Divide ice cores during these periods, and that we do not account for possible temporal variability in the correction slopes in these comparisons. The other contributing factor could be the unquantified uncertainties in the loss corrections applied to the WDC-06A and Taylor Dome records.
10. Temperature and pH effects on gas solubility and COS hydrolysis
Solubility of gases in the ocean are temperature dependent. Sea surface temperature rises during the deglaciation and this effect alone should result in an increase in gas flux out of the ocean due to decreased solubility. A + 1oC change in temperature results in less than 5% decrease in the solubilities of COS, CS2, and DMS over a range of 0 to 30 oC69–70. Present-day ocean CS2 flux correlates with temperature and can increase by about 10% per oC43; this effect is too high to be solely due to the temperature effect on CS2 solubility. Direct COS emissions are buffered against temperature increases because as solubility decreases, the hydrolysis loss rate increases at a faster rate, with + 1oC change in temperature resulting in about a 10% rise in the COS hydrolysis rate constant. The hydrolysis of COS is also sensitive to pH and slows down by about 5% per 0.1 unit decrease at pH 83,7,70. For the direct COS flux, the net impact of a + 1oC temperature increase coupled with a 0.1 unit decrease in pH is close to zero. It is not possible to drive large changes in ocean-atmosphere gas fluxes via temperature effects on gas solubility, or via the coupled effects of temperature and pH on the COS hydrolysis.
11. Calculation of the COS radiative impact based on ref. #18
The radiative impact of increased COS emissions has been estimated in chemistry-climate model experiments with two distinct emission geometries18. In the first model experiment, COS was injected into the atmosphere using the geographic distribution of present-day anthropogenic emissions at the surface, thus the emissions occur primarily over land in the northern hemisphere where the terrestrial uptake is also strongest. In the second model experiment, COS was injected at the tropical tropopause, which enables quick dispersion into the stratosphere. The location of the emissions impacts the relative increase in the tropospheric versus the stratospheric mixing ratios of COS. In the troposphere, COS acts primarily as a greenhouse gas and has a positive radiative impact. In the stratosphere, COS acts primarily as a source of sulphate aerosol with a negative radiative impact. The model also incorporates indirect radiative effects via ozone, methane, and stratospheric water vapor18.
The net (tropospheric + stratospheric) radiative impact of a 40 TgS y− 1 emission increase in the first experiment is -1.3 Wm− 2. The net radiative impact of a 6 TgS y− 1 emission increase in the second experiment is -1.5 Wm− 2. Assuming linear scaling, these sensitivities imply − 0.01 Wm− 2 (first experiment) and − 0.08 Wm− 2 (second experiment) net radiative impacts for an emission increase of 300 GgS y− 1. The deglacial increase in COS emissions is largely oceanic in origin, hence primarily in the boundary layer and in the southern hemisphere. The estimates from the two model experiments should be considered unrealistic lower and upper bounds. The best estimate may be closer to the lower bound estimate from the first model experiment.
References for Methods
60. Aydin, M., Williams, M. B., Saltzman, E. S., Feasibility of reconstructing paleoatmospheric records of selected alkanes, methyl halides, and sulfur gases from Greenland ice cores. J. Geophys. Res. 112, D07312 (2007).
61. Nicewonger, M. R., Aydin, M., Prather, M. J., Saltzman, E. S., Reconstruction of paleofire emissions over the past millennium from measurements of ice core acetylene. Geophys. Res. Lett. 47, (2020).
62. Carpenter, B. et al., Stan: A Probabilistic Programming Language. J. Statistical Software 76, 1 (2017).
63. Bowen, H. J. M., Environmental chemistry of the elements, Academic Press, (1979).
64. Fischer, H. et al., Reconstruction of millennial changes in dust emission, transport and regional sea ice coverage using the deep EPICA ice cores from the Atlantic and Indian Ocean sector of Antarctica, Earth and Planetary Science Letters, 260(1), 340–354 (2007).
65. Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., Burkner, P.-C., Rank-Normalization, Folding, and Localization: An improved \(\widehat{R}\) for assessing convergence of MCMC (with discussion). Bayesian Analysis 16, No. 2, 667–718 (2021).
66. Steig, E. et al., Wisconsinan and Holocene climate history from an ice core at Taylor Dome, Western Ross Embayment, Antarctica. Georg. Ann. 82A, 213–235 (2000).
67. Monnin, E. et al., Evidence for substantial accumulation rate variability in Antarctica during the Holocene, through synchronization of CO2 in the Taylor Dome, Dome C and DML ice cores. Earth and Plan. Sci. Lett. 224, 45–54 (2004).
68. McConnell, J. R. et al., Synchronous volcanic eruptions and abrupt climate change ~ 17.7 ka plausibly linked by stratospheric ozone depletion. Proc. Natl. Acad. Sci. 114, 10035–10040 (2017).
69. De Bruyn, W. L., Swartz, E., Hu, J. H., Shorter, J. A., Davidovits, P., Worsnop, D. R., Zahniser, M. S. and Kolb, C. E., Henry's law solubilities and Setchenow coefficients for biogenic reduced sulfur species obtained from gas-liquid uptake measurements. J. Geophys. Res. 100, 7245–7251, (1995).
70. Elliot, S., Lu, E., and Rowland, F. S., Rates and mechanisms for the hydrolysis of carbonyl sulfide in natural waters. Environ. Sci. Technol. 23, 458–461, (1989).