# Isotopic constraint on the twentieth-century increase in tropospheric ozone

Laurence Yeung, Lee. Murray, Patricia Martinerie, Emmanuel Witrant, Huanting Hu, Asmita Banerjee, Anais Orsi, Jérôme Chappellaz

► To cite this version:

Laurence Yeung, Lee. Murray, Patricia Martinerie, Emmanuel Witrant, Huanting Hu, et al.. Isotopic constraint on the twentieth-century increase in tropospheric ozone. Nature, 2019, 570 (7760), pp.224-227. 10.1038/s41586-019-1277-1 . hal-02179376

HAL Id: hal-02179376

https://hal.science/hal-02179376v1

Submitted on 8 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Isotopic constraint on the twentieth-century increase in tropospheric ozone

Laurence Y. Yeung<fnref n="1" />*, Lee. T. Murray², Patricia Martinerie<fnref n="3" />, Emmanuel Witrant⁴, Huunting Hu<fnref n="1" />,⁶, Asmita Banerjee<fnref n="1" />, Anaïs Orsi⁵ & Jérôme Chappellaz<fnref n="3" />

Tropospheric ozone ($O_3$) is a key component of air pollution and an important anthropogenic greenhouse gas<fnref n="1" />. During the twentieth century, the proliferation of the internal combustion engine, rapid industrialization and land-use change led to a global-scale increase in $O_3$ concentrations²,<fnref n="3" />; however, the magnitude of this increase is uncertain. Atmospheric chemistry models typically predict⁴⁻⁷ an increase in the tropospheric $O_3$ burden of between 25 and 50 per cent since 1900, whereas direct measurements made in the late nineteenth century indicate that surface $O_3$ mixing ratios increased by up to 300 per cent⁸⁻¹⁰ over that time period. However, the accuracy and diagnostic power of these measurements remains controversial². Here we use a record of the clumped-isotope composition of molecular oxygen ($^{18}O^{18}O$ in $O_2$) trapped in polar firn and ice from 1590 to 2016 AD, as well as atmospheric chemistry model simulations, to constrain changes in tropospheric $O_3$ concentrations. We find that during the second half of the twentieth century, the proportion of $^{18}O^{18}O$ in $O_2$ decreased by $0.03 \pm 0.02$ parts per thousand (95 per cent confidence interval) below its 1590–1958 AD mean, which implies that tropospheric $O_3$ increased by less than 40 per cent during that time. These results corroborate model predictions of global-scale increases in surface pollution and vegetative stress caused by increasing anthropogenic emissions of $O_3$ precursors⁴,⁵,¹¹. We also estimate that the radiative forcing of tropospheric $O_3$ since 1850 AD is probably less than +0.4 watts per square metre, consistent with results from recent climate modelling studies¹².

Tropospheric $O_3$ is produced by the photochemical oxidation of reduced carbon species such as methane ($CH_4$) in the presence of reactive nitrogen oxide ($NO_x$) radicals, and destroyed primarily by gas-phase reactions and deposition to Earth's surface. Changes in either $NO_x$ or reduced carbon concentrations linearly explain much of the variability in simulated tropospheric $O_3$ burdens, despite the nonlinearity in its chemistry⁴,⁵. However, the magnitude of specific terms within the $O_3$ budget are still debated, especially in ancient atmospheres, despite broad agreement in the present-day tropospheric $O_3$ burden (that is, the steady-state mass of $O_3$ in the troposphere)⁵. Moreover, concentrations of $O_3$ precursors in the preindustrial atmosphere are not known well (with the exception of $CH_4$), resulting in poor constraints on past tropospheric $O_3$. The historical latitudinal distribution of $O_3$ precursors is also poorly known¹³.

Indirect proxies are available for the past concentrations of many precursors, but these proxies are often end-products of photochemical processing (for example, nitrates for $NO_x$) with model-dependent interpretations. Other potential approaches to constraining $O_3$ precursors include measuring residual trace-gas concentrations¹¹,¹⁴ or rare-isotope abundances in particles deposited on glacial ice. These tracers, however, are generally regional in scope. In addition, many undergo post-depositional alteration and overprinting from

stratospheric processing¹⁵,¹⁶. These uncertainties limit confidence in model forecasts and hindcasts of tropospheric $O_3$.

To overcome these limitations and to constrain changes in tropospheric $O_3$ from preindustrial times to the present, we obtained a record of $^{18}O^{18}O$ in $O_2$ from polar firn and ice. This record—characterized by temporal changes in the $\Delta_{36}$ value of $O_2$ (see Methods)¹⁷—is uniquely associated with atmospheric $O_3$ photochemistry: ozone photochemistry in cold environments tends to increase $\Delta_{36}$, whereas in warmer environments it tends to decrease $\Delta_{36}$, through equilibration of the isotopic distribution in $O_2$. The tropospheric $^{18}O^{18}O$ budget is therefore described by a balance between warm (tropospheric) and cold (lower-stratospheric) contributions to $\Delta_{36}$. Ultimately, this balance allows $^{18}O^{18}O$ in tropospheric $O_2$ to respond to the tropospheric $O_3$ burden. More tropospheric $O_3$ drives $\Delta_{36}$ downwards because, when globally averaged, greater isotopic equilibration occurs at warmer temperatures¹⁷. Less tropospheric $O_3$ drives $\Delta_{36}$ upwards, towards lower-stratospheric values. Other biogeochemical factors are negligible in the sub-decadal response time of the tracer¹⁸. The change in tropospheric $\Delta_{36}$ during the twentieth century may thus distinguish between end-member scenarios for the history of tropospheric $O_3$. A two-box model of the stratosphere and troposphere that successfully explains the present-day $\Delta_{36}$ value predicts that a 41% increase in $O_3$ burden (256 Tg $O_3$ increasing to 362 Tg $O_3$) will yield a $0.05\%$ decrease in the tropospheric $\Delta_{36}$ value¹⁸. Keeping all other parameters the same, a 300% increase in $O_3$ burden would yield a decrease of about $0.1\%$ in tropospheric $\Delta_{36}$.

Air trapped in the top 105 m of polar firn collected from East Antarctica (LOCK-IN site; 3,209 m above sea level (a.s.l.)), with mean ages ranging between 1992 and 2016 AD; see Methods) has a mean $\Delta_{36}$ value of $2.01\% \pm 0.04\%$ (all uncertainties given here are $1\sigma$, $n = 38$). Ice-core $O_2$ extracted from the Greenland Ice Sheet Project 2 (GISP2-D; 3,216 m a.s.l., $n = 23$), West Antarctic Ice Sheet Divide (WDC06A; 1,759 m a.s.l., $n = 7$) and Siple Dome-94 (SDM94; 621 m a.s.l., $n = 21$) ice cores shows higher $\Delta_{36}$ values on average for mean gas ages ranging from 1590 AD to 1958 AD (Fig. 1, Extended Data Fig. 7)¹⁹⁻²¹. The mean $\Delta_{36}$ values at the three ice-coring sites are indistinguishable from each other ($P = 0.56$ or larger) and yield a cumulative mean value of $2.04\% \pm 0.03\%$ ($1\sigma$, $n = 51$). The deepest firn sample, at 107.6 m, has $\Delta_{36}$ values similar to those in ice-core gases of similar age ($\Delta_{36} = 2.08\% \pm 0.04\%$; $1\sigma$, $n = 5$), but potential contamination with younger, low-$\Delta_{36}$ air during sampling precludes its use in quantitative comparisons (see Methods). Bootstrap resampling of the data ($10^6$ subsamples) indicates that the ice-core data and the firn data above 105 m represent two different populations ($P < 0.002$), with a difference in mean $\Delta_{36}$ of $0.03\% \pm 0.02\%$ (95% confidence interval (CI); Fig. 1).

The firn data and air sampled at Cape Grim (Australia), dating back to 1978 AD¹⁸, both have mean $\Delta_{36}$ values indistinguishable from that measured in Houston in 2016–2018 AD ($P > 0.05$; Fig. 1). Using the modern mean tropospheric $\Delta_{36}$ value from Houston—$1.99\% \pm 0.02\%$

<footnote n="1">¹Department of Earth, Environmental and Planetary Sciences, Rice University, Houston, TX, USA. ²Department of Earth and Environmental Sciences, University of Rochester, Rochester, NY, USA.</footnote>

<footnote n="3">³Université Grenoble Alpes, CNRS, Institut des Géosciences de l'Environnement, Grenoble, France. ⁴Université Grenoble Alpes, CNRS, Grenoble Image Parole Signal Automatique (GIPSA-Lab), Grenoble, France. ⁵Laboratoire des Sciences du Climat et de l'Environnement, Gif-sur-Yvette, France. ⁶Present address: School of Oceanology, Shanghai Jiao Tong University, Shanghai, China.</footnote>

<footnote>*e-mail: lyeung@rice.edu</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1 |** Mean atmospheric Δ<sub>36</sub> values from different archives.</figcaption></figure>

Box-and-whisker plots (top) and kernel-smoothed bootstrap probability density distributions (bottom) of mean atmospheric Δ<sub>36</sub>. Shown values of *n* correspond to the number of measurements in the dataset used for bootstrap resampling. In the box-and-whisker plots, the vertical line inside the box is the median, the box limits are the upper and lower quartiles and the whiskers are the 10th and 90th percentiles.

(95% CI, *n* = 107)—yields a cumulative tropospheric Δ<sub>36</sub> decrease of 0.04‰ ± 0.02‰ (95% CI) in the late twentieth century. However, local effects on the firn and ice-core record—that is, potential altitudinal Δ<sub>36</sub> gradients<sup>18</sup>, synoptic transport of high-Δ<sub>36</sub> air and seasonal O<sub>3</sub> production near the ice surface<sup>22</sup>—may elevate Δ<sub>36</sub> at polar sites. Agreement across three ice-core sampling sites, however, suggests that site-specific Δ<sub>36</sub> variations are negligible. We therefore use the mean ice-to-firn Δ<sub>36</sub> decrease of 0.03‰ ± 0.02‰ in the following discussion; in either case, the measurements constrain the Δ<sub>36</sub> decrease since 1850 AD to ≤0.06‰.

To obtain a temporally resolved tropospheric history of Δ₃₆ since 1850 AD, an inversion of the firn-air and shallow ice-core Δ₃₆ data was performed using a model of gas transport in firn and trapping in ice<sup>23,24</sup> (Fig. 2a). The derived atmospheric history shows that the decrease in tropospheric Δ₃₆ begins after 1950 AD and reaches a modern-like atmospheric Δ₃₆ value after 1980 AD.

The small change in atmospheric Δ₃₆ implies an increase in the tropospheric O₃ burden of less than 40% during the twentieth century, on the basis of earlier box-modelling results¹⁸; however, a more exhaustive evaluation of the Δ₃₆ budget is necessary. A decrease in tropospheric Δ₃₆ in the twentieth century can be caused by any combination of: (i) increasing global mean tropospheric and/or stratospheric temperatures, (ii) decreased stratosphere-to-troposphere transport (STT) of air, (iii) increased rates of tropospheric O₃ photochemistry, or (iv) increased O₃ photochemistry occurring in warmer regions. We consider the impacts of these mechanisms below.

Model-mean global temperatures from the Climate Model Intercomparison Project-5 imply a temperature increase of about 1.0 °C throughout the troposphere since 1850 AD²⁵. Tropospheric temperature alone can thus only account for, at most, 0.01‰ of the observed tropospheric Δ₃₆ decrease, all else being equal. However, the stratosphere also cooled by up to 2°C, tempering the signal from tropospheric warming²⁵. The role of air STT is not easily quantified because no records exist from the early twentieth century. However, its fluxes are expected to have either stayed constant or increased—rather than decreased—during the twentieth century. Although modern reanalyses and observational constraints disagree on long-term STT trends²⁶, general circulation models robustly predict that the Brewer-Dobson circulation accelerates in response to anthropogenic activity²⁷.

<figure><img src="image_3.png" /><figcaption>**Fig. 2 |** Measurement-model comparisons for increase in tropospheric O₃ since 1850 AD. **a,** Comparison of the best-fit tropospheric history of Δ₃₆ values (shaded area, representing the 2σ-equivalent uncertainty envelope of the inverse firn/ice model results) with measurements (black and red points). Error bars represent 1 s.e.m., computed using the pooled 1σ uncertainty and the number of replicates for each sample. The firn data point at 1905 AD was excluded from the best-fit calculations because of potential contamination (see main text), but it is shown here to illustrate the consistency between the firn and ice-core data. Also shown are tropospheric Δ₃₆ predictions for 1850 AD and 2005 AD, derived from the GEOS-Chem/MERRA2 model (cyan stars, with 1σ error bars corresponding to the uncertainty in the stratospheric Δ₃₆ input; see Methods) and the mean temporal evolution of the tropospheric Δ₃₆ value derived from the GISS-E2.1 model (blue line). **b,** Hypothesis tests for the global simulations, showing compatibility with the firn and ice-core data only when the O₃ increase since 1850 AD is limited (for example, +25% for GISS-E2.1). Time-traces are five-year moving averages and shaded areas represent 1σ uncertainty in the stratospheric Δ₃₆ input. The ‘unscaled’ time trace corresponds to the GISS-E2.1 base case described in the main text. The grey box shows the Δ₃₆ change that is compatible with the firn and ice-core data. **c,** Modelled increase in tropospheric O₃ burden from the GEOS-Chem/MERRA2 (cyan stars) and GISS-E2.1 models (blue line), plotted with estimates of historical NOₓ emissions from ref. ³ (orange line).</figcaption></figure>

In particular, an increase of 1.8% per decade in the shallowest branch of the circulation has been predicted²⁸ for the time period 1960–2000, resulting in a 7% increase in STT mass flux and a 0.01‰ increase in tropospheric Δ₃₆, instead of a decrease¹⁸. The twentieth-century change in STT mass flux may therefore nullify the Δ₃₆ decrease caused by tropospheric warming.

The Δ₃₆ record thus implies that an increase in the tropospheric O₃ burden and/or the amount of O₃ photochemistry occurring in warmer tropospheric regions took place in the late twentieth century.
<figure><img src="image_4.png" /><figcaption>**Fig. 3** | Zonal-mean change in O₃ concentrations between 1850 and 2005 AD. **a,** Results obtained from simulations with the GEOS-Chem/MERRA2 (**a**) and GISS-E2.1 (**b**) models. The GEOS-Chem/MERRA2 model shows smaller O₃ concentration changes than GISS-E2.1 in the middle and upper troposphere.</figcaption></figure>

To quantify this tropospheric O₃ increase, three-dimensional global simulations were performed using the GEOS-Chem chemical-transport model driven by Modern Era Retrospective Reanalysis-2 (MERRA2) meteorology²⁹ and the Goddard Institute for Space Studies-ModelE2.1 (GISS-E2.1) chemistry-climate model³⁰. The simulation outputs provide boundary conditions for the two-box model of the stratosphere and troposphere described previously¹⁸ that translates those results to Δ₃₆ predictions.

The GEOS-Chem/MERRA2 model yields a tropospheric Δ₃₆ decrease of 0.04‰–0.05‰ between 1850 AD and 2005 AD, consistent with observations (Fig. 2a, b). This modelled Δ₃₆ decrease is associated with a 32% increase in the tropospheric O₃ burden (273 Tg O₃ increasing to 359 Tg O₃), driven primarily by increases in surface emissions of NOₓ and reduced carbon compounds⁵,³⁰ (Figs. 2c and 3). By contrast, the GISS-E2.1 model predicts a tropospheric Δ₃₆ decrease of 0.06‰–0.07‰ (computed as the mean of five ensemble members), which is generally consistent with, but slightly larger than, the observations. As implemented, the total tropospheric O₃ burden increase in the GISS-E2.1 model is 59% (221 Tg O₃ increasing to 352 Tg O₃). Some discrepancies in model-predicted Δ₃₆ values may arise from subtle inaccuracies in the strength of tropospheric convection and/or the terms of the chemical O₃ budget, which primarily affect the distribution

of O₃ in the troposphere rather than the global burden¹⁸. The larger change in tropospheric O₃ burden predicted by the GISS-E2.1 model may be related to its chemical scheme, which lacks the comprehensive halogen chemistry that is present in GEOS-Chem³¹. However, the GISS-E2.1 model reproduces the timing of the decrease in Δ₃₆ values observed in the firn and ice-core records (Fig. 2).

The GISS-E2.1 model outputs were therefore used to test hypotheses about the evolution of tropospheric O₃ and Δ₃₆ values. First, we examined whether the tropospheric Δ₃₆ decrease could be caused by changes in STT mass flux by modelling the fate of ⁷Be (half-life of 53 days) produced in the stratosphere⁴; no secular shifts in the simulated stratospheric air fraction were observed at the study sites (Extended Data Fig. 8). Second, we tested whether a 200% or 300% increase in the tropospheric O₃ burden since 1850 AD could be consistent with the data by scaling the temporal evolution of modelled tropospheric O₂ equilibration rates (Eₜrop) and keeping all other box-model inputs unchanged. The predicted decrease in tropospheric Δ₃₆ in those scenarios is outside the range consistent with observations (Fig. 2b). Third, we determined the influence of the spatial distribution of O₃ photochemistry by isolating the contributions from Eₜrop and the effective O₂ equilibration temperatures (Tₐff, with higher Tₐff reflecting more photochemistry near the surface and/or in the tropics). We found that the increase in Tₐff over the twentieth century accounts for about 60% of the decrease in tropospheric Δ₃₆, whereas the modelled evolution of Eₜrop drives about 40% of the decrease in tropospheric Δ₃₆. Therefore, the observed decrease in atmospheric Δ₃₆ was driven by increased concentrations of O₃ in the lower troposphere and towards the Equator (Fig. 3).

We conclude that the total tropospheric O₃ burden probably increased by less than 40% between 1850 AD and 2005 AD, primarily near the surface, with the majority of the increase occurring between 1950 AD and 1980 AD. A synchronous increase in oceanic emissions of iodine—the product of O₃–iodide reactions at the air–sea interface—seems to corroborate this timing, as well as the importance of halogen chemistry in late twentieth-century O₃ budgets³¹,³². Moreover, historical global emission estimates of O₃ precursors³, combined with current model chemical schemes, appear to capture the main features of the tropospheric O₃ increase since 1850 AD. We stress that these findings do not invalidate the late nineteenth-century surface O₃ observations, which may simply not be representative of the tropospheric O₃ burden at that time.

These results also constrain the anthropogenic influence on the radiative forcing of tropospheric O₃. The GEOS-Chem/MERRA2 model yields an instantaneous radiative forcing at the tropopause due to tropospheric O₃ of +0.33 W m⁻² between 1850 and 2005 AD, consistent with the multi-model mean calculated for the most recent iteration of the Climate Model Intercomparison Project (CMIP6; +0.31 ± 0.06 W m⁻²)¹². The GISS-E2.1 model yields a slightly higher tropospheric O₃ radiative forcing of +0.42 W m⁻² from 1850 AD to the 2000s AD (Extended Data Fig. 9) but also tends to overestimate the observed decrease in Δ₃₆. The greater ability of GEOS-Chem to represent historical tropospheric O₃ chemistry may translate to more accurate predictions of past and future tropospheric O₃ changes and radiative forcing, particularly in regions where O₃ precursor emissions are changing rapidly.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-019-1277-1.

Received: 11 September 2018; Accepted: 10 May 2019; Published online 12 June 2019.

1. National Research Council. *Rethinking the Ozone Problem in Urban and Regional Air Pollution* (National Academies Press, 1991).

2. Cooper, O. R. et al. Global distribution and trends of atmospheric ozone: an observation-based review. *Elementa* **2**, 000029 (2014).
3. Lamarque, J. F. et al. Historical (1850–2000) gridded anthropogenic and biomass burning emissions of reactive gases and aerosols: methodology and application. *Atmos. Chem. Phys.* **10**, 7017–7039 (2010).

4. Murray, L. T. et al. Factors controlling variability in the oxidative capacity of the troposphere since the Last Glacial Maximum. *Atmos. Chem. Phys.* **14**, 3589–3622 (2014).

5. Young, P. J. et al. Pre-industrial to end 21st century projections of tropospheric ozone from the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP). *Atmos. Chem. Phys.* **13**, 2063–2090 (2013).

6. Parrella, J. P. et al. Tropospheric bromine chemistry: implications for present and pre-industrial ozone and mercury. *Atmos. Chem. Phys.* **12**, 6723–6740 (2012).

7. Archibald, A. T. et al. Impacts of HOₓ regeneration and recycling in the oxidation of isoprene: Consequences for the composition of past, present, and future atmospheres. *Geophys. Res. Lett.* **38**, L05804 (2011).

8. Volz, A. & Kley, D. Evaluation of the Montsouris series of ozone measurements made in the nineteenth century. *Nature* **332**, 240–242 (1988).

9. Marenco, A., Gouget, H., Nédélec, P. & Pagès, J.-P. Evidence of a long-term increase in tropospheric ozone from Pic du Midi data series: consequences: positive radiative forcing. *J. Geophys. Res.* **99**, 16617–16632 (1994).

10. Pavelin, E. G., Johnson, C. E., Rughoopth, S. & Toumi, R. Evaluation of pre-industrial surface ozone measurements made using Schönbein's method. *Atmos. Environ.* **33**, 919–929 (1999).

11. Worton, D. R. et al. Evidence from firn air for recent decreases in non-methane hydrocarbons and a 20th century increase in nitrogen oxides in the northern hemisphere. *Atmos. Environ.* **54**, 592–602 (2012).

12. Checa-Garcia, R., Hegglin, M. I., Kinnison, D., Plummer, D. A. & Shine, K. P. Historical tropospheric and stratospheric ozone radiative forcing using the CMIP6 database. *Geophys. Res. Lett.* **45**, 3264–3273 (2018).

13. Zhang, Y. et al. Tropospheric ozone change from 1980 to 2010 dominated by equatorward redistribution of emissions. *Nat. Geosci.* **9**, 875–879 (2016).

14. Newland, M. J. et al. Changes to the chemical state of the Northern Hemisphere atmosphere during the second half of the twentieth century. *Atmos. Chem. Phys.* **17**, 8269–8283 (2017).

15. Alexander, B. et al. Quantifying atmospheric nitrate formation pathways based on a global model of the oxygen isotope composition (Δ¹⁷O) of atmospheric nitrate. *Atmos. Chem. Phys.* **9**, 5043–5056 (2009).

16. Erbland, J. et al. Air-snow transfer of nitrate on the East Antarctic Plateau – part 1: isotopic evidence for a photolytically driven dynamic equilibrium in summer. *Atmos. Chem. Phys.* **13**, 6403–6419 (2013).

17. Yeung, L. Y., Ash, J. L. & Young, E. D. Rapid photochemical equilibration of isotope bond ordering in O₂. *J. Geophys. Res.* **119**, 10552–10566 (2014).

18. Yeung, L. Y. et al. Isotopic ordering in atmospheric O₂ as a tracer of ozone photochemistry and the tropical atmosphere. *J. Geophys. Res. Atmos.* **121**, 12541–12559 (2016).

19. Taylor, K. C. et al. Dating the Siple Dome (Antarctica) ice core by manual and computer interpretation of annual layering. *J. Glaciol.* **50**, 453–461 (2004).

20. Vinther, B. M. et al. A synchronized dating of three Greenland ice cores throughout the Holocene. *J. Geophys. Res. Atmos.* **111**, (2006).

21. Buizert, C. et al. The WAIS Divide deep ice core WD2014 chronology – part 1: methane synchronization (68–31 ka BP) and the gas age–ice age difference. *Clim. Past* **11**, 153–173 (2015).

22. Legrand, M. et al. Inter-annual variability of surface ozone at coastal (Dumont d'Urville, 2004–2014) and inland (Concordia, 2007–2014) sites in East Antarctica. *Atmos. Chem. Phys.* **16**, 8053–8069 (2016).

23. Witrant, E. & Martinerie, P. Input estimation from sparse measurements in LPV systems and isotopic ratios in polar firns. *IFAC Proc.* **46**, 659–664 (2013).

24. Witrant, E. et al. A new multi-gas constrained model of trace gas non-homogeneous transport in firn: evaluation and behaviour at eleven polar sites. *Atmos. Chem. Phys.* **12**, 11465–11483 (2012).

25. Santer, B. D. et al. Human and natural influences on the changing thermal structure of the atmosphere. *Proc. Natl Acad. Sci. USA* **110**, 17235–17240 (2013).

26. Boothe, A. C. & Homeyer, C. R. Global large-scale stratosphere-troposphere exchange in modern reanalyses. *Atmos. Chem. Phys.* **17**, 5537–5559 (2017).

27. Butchart, N. et al. Chemistry-climate model simulations of twenty-first century stratospheric climate and circulation changes. *J. Clim.* **23**, 5349–5374 (2010).

28. Lin, P. & Fu, Q. Changes in various branches of the Brewer–Dobson circulation from an ensemble of chemistry climate models. *J. Geophys. Res. Atmos.* **118**, 73–84 (2013).

29. Gelaro, R. et al. The Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2). *J. Clim.* **30**, 5419–5454 (2017).

30. Schmidt, G. A. et al. Configuration and assessment of the GISS ModelE2 contributions to the CMIP5 archive. *J. Adv. Model. Earth Syst.* **6**, 141–184 (2014).

31. Sherwen, T., Evans, M. J., Carpenter, L. J., Schmidt, J. A. & Mickley, L. J. Halogen chemistry reduces tropospheric O₃ radiative forcing. *Atmos. Chem. Phys.* **17**, 1557–1569 (2017).

32. Legrand, M. et al. Alpine ice evidence of a three-fold increase in atmospheric iodine deposition since 1950 in Europe due to increasing oceanic emissions. *Proc. Natl Acad. Sci. USA* **115**, 12136–12141 (2018).

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2019
## METHODS

**Isotopic measurements.** The $\Delta_{36}$ tracer characterizes the proportions of $^{18}$O$^{18}$O compared to the stochastic distribution of isotopes in O$_2$, that is,

$$ \Delta_{36} = \left( \frac{^{36}R_{\text{measured}}}{^{36}R_{\text{stochastic}}} - 1 \right) \qquad (1) $$

$$ ^{36}R_{\text{measured}} = \frac{[^{18}\text{O}^{18}\text{O}]}{[^{16}\text{O}^{16}\text{O}]} \quad ^{36}R_{\text{stochastic}} = \frac{[^{18}\text{O}][^{18}\text{O}]}{[^{16}\text{O}][^{16}\text{O}]} = (^{18}\text{R})^2 \quad (2) $$

where $^{18}$R is the bulk $^{18}$O/$^{16}$O ratio in O$_2$. The $\Delta_{36}$ value, reported in parts per thousand (‰), reflects proportional deviations from the random occurrence of isotopes in O$_2$ molecules$^{33–35}$: $\Delta_{36} = 0$ represents a random distribution of isotopes for a given $^{18}$R.

Firn-air samples for $\Delta_{36}$ analysis (70–100 µmol O$_2$) were purified using a gas chromatography system described previously$^{18,36}$. Ice-core samples (70–110 g) were extracted using the wet-extraction technique$^{37}$. The exterior of the ice was trimmed by about 5 mm before placing the sample in a vessel that was attached to a turbo-pumped vacuum line. The bottom of the vessel was placed in an ethanol cold bath between –30°C and –40°C and attached to the line at a KF joint before the headspace was evacuated and turbo-pumped for 30 min. The sample was then melted in warm water while stirring (glass stir rod) while the air was condensed, after passing through two cryogenically cooled U-traps (–196°C), onto a degassed silica gel adsorbent at –196°C. Samples were then transferred to the gas chromatography line.

Firn and ice-core $\Delta_{36}$ analyses were run separately by three analysts over several years to mitigate temporal and user bias. The pooled standard deviation for replicate measurements of air-O$_2$ in Houston was ±0.031‰ and ±0.044‰ for $^{18}$O and $\Delta_{36}$, respectively (n = 107). The pooled standard deviation for replicate firn-air measurements was ±0.026‰ and ±0.039‰ for $^{18}$O and $\Delta_{36}$, respectively (14 depths). The pooled standard deviation for replicate ice-core measurements was ±0.044‰ (typically better than ±0.03‰, with several exceptions) and ±0.029‰ for $^{18}$O and $\Delta_{36}$, respectively (22 depths). Absolute $\Delta_{36}$ values were calibrated during each analytical session (that is, the lifetime of a mass spectrometer filament) using photochemical equilibration and gas heating$^{17}$. The absolute accuracy in $\Delta_{36}$ values was estimated to be ±0.02‰, limited by the temperature stability of the calibration experiments (±1 K). Isotopic measurements of the $^{61}$N value of N$_2$ and the $^{18}$O value of O$_2$ in the firn-air samples were performed at Laboratoire des Sciences du Climat et de l'Environnement, France (LSCE) using published methods$^{38}$.

**Atmospheric trend estimates from firn-air and ice-core data.** Gas diffusion through the open porosity of firn, followed by porosity closure, results in increasing gas ages with depth$^{39,40}$ and gas mixing between layers. Trace-gas records in ice cores are further mixed by progressive gas enclosure into individual bubbles or layered gas trapping$^{41}$. These processes are often taken into account in a simplified way in ice-core records by attributing a mean gas age to samples and correcting for gravitational effects. Here, we use a model of firn physics that allows a more precise reconstruction of atmospheric variations.

*Model of gas transport in firn and trapping in ice.* The IGE-GIPSA model has been described previously$^{24}$. It includes gravitational settling, molecular diffusion and advection of air together with firn sinking. 'Convection' (gas mixing due to pressure/wind effects) is considered only in the upper firn. Thermal diffusion (mostly affecting the upper firn) and molecular size-dependent fractionation$^{42,43}$ are neglected. The effect of deep firn fractionations for O$_2$ is not fully understood$^{44}$. The IGE-GIPSA model ranked second among six models in a firn-model inter-comparison$^{40}$, with nearly equal performance with the first-ranked model for model-data discrepancies.

Recent improvements$^{23}$ to the IGE-GIPSA inverse model provide a more robust definition of the optimal solution and a more rigorous treatment of isotopes. However, the new isotopic model has not yet been extended to ice; hence, a simplified isotopic model$^{45}$, based on firn-fractionation correction$^{46}$, was used for the multi-site-constrained trend reconstructions. The results of the simplified and new isotopic models are compared for LOCK-IN firn below.

*Model tuning to the drilling sites used.* Firn-air sampling at the LOCK-IN site followed the principles first described in Schwander et al.$^{39}$. The site, located 136 km away from the Concordia Antarctic station towards Dumont D'Urville (74° 8.310' S, 126° 9.510' E; 3,209 m a.s.l.), was documented in January 2016 initially to investigate firn physics$^{47}$. Its mean annual temperature (borehole temperature at a depth of 20 m) is –53.15°C and the best-estimate accumulation rate is 3.6 cm water equivalent per year. A 3-m-long Tec-Flex 40 rubber bladder was used to seal the borehole at each sampled depth. A 'Bender baffle' was attached to the bladder for venting firn-air in direct contact with the bladder, while the sample line below the baffle was directed towards gas analysers and canisters. 0.25-inch

Dekabon tubing was used between the bladder and the surface. CO$_2$, CO and CH$_4$ concentrations were measured in the field (using LiCor and SARA optical analysers, calibrated daily against Scott-Marrin gas standards) in Silcoan canisters and in CSIRO glass flasks, allowing us to check for contamination. Consistent mixing ratios were obtained, except at the greatest depth where $\Delta_{36}$ was measured (107.65 m), which could not be sampled correctly, probably owing to contamination amplified by the reduced flow rate from deep firn (about 0.04 L min$^{-1}$). Firn diffusivity was constrained using CO$_2$ and CH$_4$ data obtained either in the field or in canisters at Institut des Géosciences de l'Environnement, Grenoble (IGE), and by CFC-11, CFC-12, CFC-113, HFC-134a and SF$_6$ data measured at the University of East Anglia. We note that because physical fractionation between $^{16}$O$^{18}$O and $^{18}$O$^{18}$O depends on mass, the relative proportions of the two species, and therefore the $\Delta_{36}$ values, have a weak sensitivity to physical fractionation mechanisms in the firn$^{48}$.

The model was tuned to the Summit site in Greenland, where the GISP2-D ice core was obtained and to Siple Dome$^{24}$. However, the Siple Dome firn-air pumping site$^{49}$ has a different accumulation rate from the SDM94 drilling site (from which samples for $\Delta_{36}$ measurements were obtained)$^{19,50}$. Owing to the weak sensitivity of $\Delta_{36}$ to firn fractionations$^{48}$, this disparity probably has little effect on the results; however, test simulations with and without Siple data were performed. The IGE-GIPSA model has not been tuned to the West Antarctic Ice Sheet (WAIS) Divide. However, the North Greenland Eemian Ice Drilling (NEEM) site has a similar temperature and accumulation rate$^{24,40}$, so the NEEM firn characteristics were used for the WAIS Divide.

The firn physics model works with a constant accumulation rate so the depth-age relationship cannot exactly follow the ice-core chronology. To preserve the correct gas chronology, the ice sample depths were shifted so the best-estimate mean gas age of the samples matched the model-mean gas age. The resulting depth shifts are small for GISP2 (0 to 1.24 m) and the WAIS Divide (1.3 to 2 m) but larger for Siple Dome (5 to 5.8 m). The modelled O$_2$ age distributions (Extended Data Fig. 1) show that $^{16}$O$^{16}$O and $^{18}$O$^{18}$O ages differ by less than a year. Age distributions in ice are primarily controlled by accumulation rates: they control the speed at which a firn layer sinks through the porosity closure zone$^{41}$. However, other firn physics parameters also influence the gas age distributions. For example, the blue lines in Extended Data Fig. 1 illustrate the effect of using firn diffusivity tuned to the NEEM-EU or NEEM-US air pumping sites$^{24,40}$. Age distributions in ice are more uncertain than in firn because the porosity closure process is weakly constrained (for example, ref. 47). Our age distributions are roughly consistent with published results: somewhat wider than those of Mitchell et al.$^{51}$ for the WAIS Divide but less smooth than those of Ahn et al.$^{52}$ for Siple Dome.

Extended Data Fig. 2 compares the results of our most precise firn model (written in isotopic $\delta$ notation$^{23}$) and the simplified model used to include ice data$^{45}$. All results are within uncertainty limits of one another. The preferred solutions, obtained without constraining the models with the contaminated sample (107.65 m), are nearly identical. The solutions obtained when using the 107.65-m sample appear more divergent because the two models select a slightly different optimal solution. This difference is insignificant. We note that the choice of the optimal solution is an important issue for our inverse problem, which has an infinite number of solutions$^{53}$. A statistical robustness-oriented optimal solution, adapted to small datasets, was implemented$^{23}$.

*$\Delta_{36}$ for firn modelling.* The species-specific characteristics used in the firn model are the molecular mass and diffusion coefficient, D. The relative diffusion coefficients with respect to CO$_2$ for O$_2$ isotopologues were calculated as in the supplement of Buizert et al.$^{40}$: $D(O_2)/D(CO_2) = 1.302$, $D(^{16}\text{O}^{18}\text{O})/D(CO_2) = 1.2837$ and $D(^{18}\text{O}^{18}\text{O})/D(CO_2) = 1.2668$. Using $\Delta_{36}$ makes the value dependent on three species ($^{18}\text{O}^{18}\text{O}$, $^{16}\text{O}^{18}\text{O}$ and $^{16}\text{O}^{16}\text{O}$) rather than two with $\delta$ notation ($^{18}\text{O}^{18}\text{O}$ and $^{16}\text{O}^{16}\text{O}$). Although $\Delta$ notation requires more information, it is insensitive to fractionation processes in firn because mass-dependent fractionations, such as gravitational fractionation, have minimal effects on $\Delta_{36}$ values$^{48}$. Extended Data Fig. 3 illustrates gravitational fractionation affecting the $^{61}$N value of N$_2$ and the $^{18}$O and $^{18}\text{O}^{18}\text{O}$ values of O$_2$, while leaving $\Delta_{36}$ unaffected. The deviations of $^{61}$N, $^{18}\text{O}^{18}\text{O}$ and $^{18}\text{O}^{18}\text{O}/4$ from the barometric slope are correlated. Even for the anomalously low values of the sample taken at a depth of 94 m, for which atmospheric contamination is likely, the correlation between $\delta$ values remains, and $\Delta_{36}$ shows no anomaly (Extended Data Fig. 4).

The firn model cannot directly work with $\Delta$ values without major modifications. Therefore we converted $\Delta_{36}$ into the gravitationally corrected $^{61}\text{O}^{18}\text{O}$ and ran the firn model with a null gravity constant (after Buizert et al.$^{40}$). $\Delta_{36}$ values are converted to $^{61}\text{O}^{18}\text{O}$ values using

$$ ^{61}\text{O}^{18}\text{O} = (\Delta_{36} + 1)(^{61}\text{O}^{18}\text{O} + 1)^2 - 1 \quad (3) $$

Assuming that $^{61}\text{O}$ is constant in the atmosphere on centennial timescales, that is, $^{61}\text{O}_{\text{atm}} = 0$, $\Delta_{36}$ is equal to $^{61}\text{O}^{18}\text{O}$ in the absence of firn fractionation.
Extended Data Fig. 4 shows that Δ36 and gravity-corrected δ¹⁸O¹⁸O values are similar. However, the method used for gravity correction is important: because of correlated analytical biases (Extended Data Fig. 3), using δ¹⁸O measurements made at Rice University (δ¹⁸O<sub>Rice</sub>) for both the Δ36-to-δ¹⁸O¹⁸O conversion and the gravity correction (δ¹⁸O¹⁸O<sub>grav.cor.</sub> = δ¹⁸O¹⁸O - 2 × δ¹⁸O<sub>Rice</sub>) provides the most consistent values with the initial Δ36 signal (Extended Data Fig. 4). Using δ¹⁸O data to perform the gravity correction implicitly assumes that deep-firn fractionation is negligible, or that the deep-firn fractionation mechanism preserves Δ36 values. No sign of systematic fractionation is observed in the δ¹⁸O LOCK-IN data, although the higher variability in the 90–110 m range (Extended Data Fig. 3) is unexplained. Nevertheless, for plausible physical isotope fractionation mechanisms stemming from molecular diffusion, the larger variability in δ¹⁸O (typically 0.1‰, but as large as 0.3‰ at 94 m) would result in Δ36 errors not exceeding 0.015‰. Δ36 values are thus probably unaffected by fractionation at porosity closure.

Contamination with modern atmospheric O₂ is not observable in the Δ36 values of these samples, partly because the magnitude of the secular shift is small; for example, a 30% contamination of modern air in pre-1950 air would yield a Δ36 error of 0.015‰ if they differed by 0.05‰. The agreement in Δ36 values between pristine and potentially contaminated depths corroborates the expectation that such errors are not detectable.

Large O₂ fractionations are observed on ice cores stored for many years in cold rooms, so very cold storage (about -50°C) is required to preserve the O₂/N₂ ratio⁵⁴. Our Δ36 measurements in ice were performed on ice cores drilled at very different dates—between 1990 (GISP2) and 2007 (WAIS Divide)—and show consistent results.

An ¹⁶O¹⁶O concentration trend is necessary to reconstruct the atmospheric trend in δ¹⁸O¹⁸O or Δ36 (ref. ²³). To test its impact on our results, two scenarios were used: a constant-O₂ trend and a databases trend (SCRIPPS; http://scrippso2.ucsd.edu/) extrapolated for dates before 1991 based on a constant O₂/CO₂ ratio (-11.2 per meg per p.p.m. CO₂ in the atmosphere). For the SCRIPPS-based scenario, the O₂ trend induces a diffusional fractionation in firn, affecting both δ¹⁸O and δ¹⁸O¹⁸O. The δ¹⁸O data corrected for this fractionation (that is, δ¹⁸O<sub>Rice</sub> - δ¹⁸O<sub>difus.</sub>)⁴⁶,⁵⁵ are then used to correct δ¹⁸O¹⁸O values for gravitational fractionation. The diffusional fractionation of δ¹⁸O¹⁸O is taken into account when using the simplified method for isotopes⁴⁵. The Δ36 trends obtained with the SCRIPPS-based or constant-O₂ trends are nearly identical (Extended Data Fig. 5); thus the atmospheric trend in O₂ does not affect Δ36.

Alternatively, O₂, δ¹⁸O and guessed Δ36 atmospheric trends can be converted into ¹⁶O¹⁶O, ¹⁶O¹⁸O and ¹⁸O¹⁸O trends for input into the forward firn model (with gravity). The results (Extended Data Fig. 5) show that using a constant atmospheric value for Δ36 yields nearly constant values in firn and ice. This confirms the weak sensitivity of Δ36 to firn physics. Using the Δ36 trends from the atmospheric box model (and a constant value before 1850) shows results consistent with Fig. 2: the +25% scenario remains within error of most ice-core data whereas the +200% and +300% scenarios overestimate the ice-core data. All scenarios are consistent with firn data.

Sensitivity tests on atmospheric trend reconstructions. The preferred trend (excluding the contaminated LOCK-IN firn depth of 107.65 m) is shown in black in Extended Data Fig. 5. Sensitivity tests were performed to evaluate the impact of uncertain parameters on the reconstructed Δ36 trend. The analysis of the LOCK-IN firn and ice is ongoing; thus uncertainties were tested by using different accumulation rates (3 to 4.5 cm yr⁻¹), density profiles (based on field volume and weight measurements or γ-ray absorption measurements performed at Alfred Wegener Institute, Bremerhaven), closed-to-open porosity parameterizations (equations (4)–(6) in Witrant et al.²⁴ with the full close-off depth z<sub>F</sub> varying between 109 and 117 m), and two (CO₂ and CH₄) or six reference gases for diffusivity calculation. As long as a firn diffusivity profile is calculated consistently with the parameters used for the trend reconstruction, they have almost no impact on the Δ36 reconstructed trend (Extended Data Fig. 6, right panel).

Simulations excluding Siple Dome or WAIS Divide data, and representing WAIS Divide as the NEEM-US site rather than the NEEM-EU site, were performed (Extended Data Fig. 6, left panel). The results remain within uncertainty limits and only the early part of the trend is affected. Multiple relative minima in the cost function are unusual with the new definition of the optimal solution²³, but their occurrence in the Δ36 simulations may be due to the scarcity of ice data and relatively high noise in the data. Systematic searches of the relative minima were performed, and the related solutions also remain close to the preferred solution except in the earliest part of the trend.

**Chemical-transport and chemistry-climate modelling.** Atmospheric simulations. Two simulations were performed using v11-02c of the GEOS-Chem global three-dimensional chemical-transport model to estimate the impact of anthropogenic emissions on the tropospheric O₃ burden. The chemistry was determined with a unified chemical mechanism from the surface to the mesopause⁵⁶. A present-day simulation was performed using meteorology from the MERRA-2 reanalysis²⁹

and all default emissions for the year 2005, which was compared to an identical simulation in which all anthropogenic combustion and fertilizer sources were removed. The simulations were initialized over the same year of meteorology for nine iterations to allow the STT flux to equilibrate, with the tenth iteration used for analysis. Diagnostics necessary for calculating Δ36 (mass; specific humidity; temperature; concentrations of O₃, O(¹D) and O(³P); and tropopause pressure)¹⁸ were archived at monthly resolution.

Five transient ensemble-member model runs using the GISS-E2.1 model with prescribed climatological sea-surface temperatures and fully coupled online reactive chemistry were performed using CMIP6 boundary conditions in 1850–2015 ad. The same diagnostics as those used in GEOS-Chem were archived at monthly resolution. Five ensemble members were averaged to account for internal climate variability. A separate ⁷Be transient run was implemented and performed with decadal short-lived forcings prescribed from the fully coupled ensemble runs. The ⁷Be source parameterization⁵⁷ reflects changes in solar activity and Earth's magnetic field.

Monthly averaged model outputs were used to derive E<sub>trop</sub> and T<sub>eff</sub> values using methods described previously¹⁸. A chemical threshold (150 parts per billion O₃) was used to define the tropopause¹⁸ at each monthly time-step. The specific concentration threshold has minimal impact on the modelled Δ36 value (tropospheric variations of about 0.01‰; see ref. ¹⁸), but representing STT fluxes is important. Using the World Meteorological Organization's lapse-rate tropopause, which may not be uniquely defined in the subtropics, results in poorer agreement between the model and Δ36 measurements (for example, Δ36 = 2.05‰ modelled for the present) and a larger tropospheric O₃ burden. The resulting calculated increase in tropospheric O₃ burden in 1850–2015 ad is smaller (+48% for GISS-E2.1; 278–412 Tg O₃).

Mean tropospheric Δ36 values were computed from model outputs by solving the family of mass-balance equations for O₂ isotopologues, for example⁴⁸

$$F_{\text{ST}}\chi_S^3 - F_{\text{TS}}\chi_T^3 + F_P\chi_P^3 - F_R\chi_T^3 \alpha_{36,R} - E_{\text{trop}}(\chi_T^3 - \chi_{\text{Tequil}}^3) = 0 \quad (4)$$

The first two terms are isotopic fluxes from the stratosphere into the troposphere (F<sub>ST</sub>χ<sub>S</sub>³ and from the troposphere into the stratosphere (F<sub>TS</sub>χ<sub>T</sub>³). The mass fluxes F<sub>ST</sub> and F<sub>TS</sub> are assumed to be equal on the timescales relevant to this work. The tropospheric isotopic mole fraction is described by χ<sub>T</sub>³. The stratospheric mixing endmember (χ<sub>S</sub>³) is Δ36,<sub>strat</sub> = 2.32‰ ± 0.08‰; the uncertainty gives rise to the range in predicted post-1950 decreases in tropospheric Δ36 documented in the main text. The next two terms, F<sub>P</sub>χ<sub>P</sub>³ and F<sub>R</sub>χ<sub>T</sub>³, which are nearly negligible on these timescales, correspond to O₂ production by photosynthesis and O₂ consumption via respiration, respectively⁵⁸. The final term, E<sub>trop</sub> (χ<sub>T</sub>³ - χ<sub>T<sub>equil</sub></sub>³), describes isotopic reordering towards Δ36,<sub>tequil</sub> (the equilibrium Δ36 value for T<sub>eff</sub>)¹⁷,³⁵. On the timescales relevant to this study, equation (4) can be approximated within ~0.01‰ by

$$F_{\text{ST}}(\Delta_{36,\text{strat}} - \Delta_{36,\text{trop}}) - E_{\text{trop}}(\Delta_{36,\text{trop}} - \Delta_{36,\text{Tequil}}) = 0 \quad (5)$$

Equation (4) was solved in monthly time-steps using the outputs of the transient simulations to update the E<sub>trop</sub>, F<sub>ST</sub>, F<sub>TS</sub> and Δ36,<sub>tequil</sub> values after a 100-year spin-up using the mean 1850 ad model inputs. The sensitivity of this model to factors other than rates of gas-phase photochemistry has been explored previously¹⁸,⁴⁸,⁵⁸. ⁷Be tests for STT. The original description of the ⁷Be tracer implementation and its application to STT in ModelE is provided in Koch et al.⁵⁹, with subsequent updates⁶⁰,⁶¹. We further updated the ⁷Be source to be a space- and time-varying function of atmospheric pressure, solar modulation potential (φ, in gigavolts) and Earth's magnetic field using the parameterization of Usoskin et al.⁵⁷.

We ran the GISS-E2.1 model with the online ⁷Be tracer in the configuration described above, but instead of including online reactive photochemistry, we prescribe monthly mean O₃ and aerosol forcings archived from the fully coupled ensemble runs. The model was initialized using 1850 ad forcings for 10 years, and then run in 1850–2016 ad, updating the forcings each month.

Historical monthly mean reconstructions of φ from ground-based neutron monitors and ionization chambers exist since July 1936 ad⁶². For earlier years, variability in φ was estimated using a linear relationship with monthly mean sunspot number SN

$$\phi = 358.6 + 2.58 \text{ SN} \quad (6)$$

using SN values provided by the Royal Observatory of Belgium (http://sidc.be/silso)⁶³. Over the period 1936–2016 ad, this relationship captures roughly half of the variability seen in the reconstruction (Extended Data Fig. 10; R² = 0.48; n = 966 months). We used the historical monthly mean magnetic North Pole reconstruction⁶⁴ provided for 1590 ad to present by NOAA through the National Centers for Environmental Information (https://www.ngdc.noaa.gov/geomag/GeomagneticPoles.shtml).
**Radiative forcing calculations.** For GISS E2.1, radiative forcing was estimated as the difference in the all-sky instantaneous radiative flux at the tropopause due to tropospheric O₃ between the 2000s and 1850s ensemble climatologies. It was determined online within the model by re-running the radiation-transfer code with tropospheric O₃ removed. The instantaneous radiative imbalance at the tropopause is a common approximation of the top-of-the-atmosphere radiative flux after allowing for fast feedbacks in formal use by international climate assessments⁶⁵,⁶⁶. For GEOS-Chem, we prescribed the monthly mean tropospheric O₃ distributions for 2005 versus 1850 AD within a present-day GISS E2.1 climate simulation using 2000s boundary conditions. The radiative forcing was again approximated as the decadal net change in the all-sky instantaneous radiative flux at the tropopause.

## Data availability

The isotopic data and main LOCK-IN firn data that support the findings of this study are available from the PANGAEA database (https://doi.pangaea.de/10.1594/PANGAEA.901154). The LOCK-IN firn analysis is ongoing, so additional firn data underlying sensitivity tests in Extended Data Fig. 6 will be published elsewhere and made available freely and immediately upon request.

## Code availability

The computer codes used to support the findings of this study are available from the authors upon reasonable request.

33. Eiler, J. M. "Clumped-isotope" geochemistry—the study of naturally-occurring, multiply-substituted isotopologues. *Earth Planet. Sci. Lett.* **262**, 309–327 (2007).

34. Eiler, J. M. & Schauble, E. ¹⁸O¹³C¹⁶O in the Earth's atmosphere. *Geochim. Cosmochim. Acta* **68**, 4767–4777 (2004).

35. Wang, Z., Schauble, E. A. & Eiler, J. M. Equilibrium thermodynamics of multiply substituted isotopologues of molecular gases. *Geochim. Cosmochim. Acta* **68**, 4779–4797 (2004).

36. Yeung, L. Y., Hayles, J. A., Hu, H., Ash, J. L. & Sun, T. Scale distortion from pressure baselines as a source of inaccuracy in triple-isotope measurements. *Rapid Commun. Mass Spectrom.* **32**:1811-1821 (2018).

37. Severinghaus, J. P., Grachev, A., Luz, B. & Caillon, N. A method for precise measurement of argon 40/36 and krypton/argon ratios in trapped air in polar ice with applications to past firn thickness and abrupt climate change in Greenland and at Siple Dome, Antarctica. *Geochim. Cosmochim. Acta* **67**, 325–343 (2003).

38. Capron, E. et al. Synchronising EDML and NorthGRIP ice cores using ⁸¹⁰O of atmospheric oxygen (⁸¹⁰O_atm) and CH₄ measurements over MIS5 (80–123 kyr). *Quat. Sci. Rev.* **29**, 222–234 (2010).

39. Schwander, J. et al. The age of the air in the firn and the ice at Summit, Greenland. *J. Geophys. Res. Atmos.* **98**, 2831–2838 (1993).

40. Buizert, C. et al. Gas transport in firn: multiple-tracer characterisation and model intercomparison for NEEM, Northern Greenland. *Atmos. Chem. Phys.* **12**, 4259–4277 (2012); corrigendum **14**, 3571–3572 (2014).

41. Fourteau, K. et al. Analytical constraints on layered gas trapping and smoothing of atmospheric variability in ice under low-accumulation conditions. *Clim. Past* **13**, 1815–1830 (2017).

42. Huber, C. et al. Evidence for molecular size dependent gas fractionation in firn air derived from noble gases, oxygen, and nitrogen measurements. *Earth Planet. Sci. Lett.* **243**, 61–73 (2006).

43. Severinghaus, J. P. & Battle, M. O. Fractionation of gases in polar ice during bubble close-off: new constraints from firn air Ne, Kr and Xe observations. *Earth Planet. Sci. Lett.* **244**, 474–500 (2006).

44. Battle, M. O. et al. Controls on the movement and composition of firn air at the West Antarctic Ice Sheet Divide. *Atmos. Chem. Phys.* **11**, 11007–11021 (2011); corrigendum **14**, 9511 (2014).

45. Wang, Z. et al. The isotopic record of Northern Hemisphere atmospheric carbon monoxide since 1950: implications for the CO budget. *Atmos. Chem. Phys.* **12**, 4365–4377 (2012).

46. Trudinger, C. M. et al. Modeling air movement and bubble trapping in firn. *J. Geophys. Res. Atmos.* **102**, 6747–6763 (1997).

47. Burr, A. et al. Pore morphology of polar firn around closure revealed by X-ray tomography. *Cryosphere* **12**, 2481–2500 (2018).

48. Yeung, L. Y., Young, E. D. & Schauble, E. A. Measurements of ¹⁸O¹⁸O and ¹⁷O¹⁸O in the atmosphere and the influence of isotope-exchange reactions. *J. Geophys. Res.* **117**, D18306 (2012).

49. Butler, J. H. et al. A record of atmospheric halocarbons during the twentieth century from polar firn air. *Nature* **399**, 749–755 (1999).

50. Kreutz, K. J., Mayewski, P. A., Whitlow, S. I. & Twickler, M. S. Limited migration of soluble ionic species in a Siple Dome, Antarctica, ice core. *Ann. Glaciol.* **27**, 371–377 (1998).

51. Mitchell, L. E. et al. Observing and modeling the influence of layering on bubble trapping in polar firn. *J. Geophys. Res. Atmos.* **120**, 2558–2574 (2015).

52. Ahn, J., Brook, E. J. & Buizert, C. Response of atmospheric CO₂ to the abrupt cooling event 8200years ago. *Geophys. Res. Lett.* **41**, 604–609 (2014).

53. Rommelaere, V., Arnaud, L. & Barnola, J.-M. Reconstructing recent atmospheric trace gas concentrations from polar firn and bubbly ice data by inverse methods. *J. Geophys. Res. Atmos.* **102**, 30069–30083 (1997).

54. Extier, T. et al. On the use of ⁸¹O_atm for ice core dating. *Quat. Sci. Rev.* **185**, 244–257 (2018).

55. Sapart, C. J. et al. Can the carbon isotopic composition of methane be reconstructed from multi-site firn air measurements? *Atmos. Chem. Phys.* **13**, 6993–7005 (2013).

56. Eastham, S. D., Weisenstein, D. K. & Barrett, S. R. H. Development and evaluation of the unified tropospheric-stratospheric chemistry extension (UCX) for the global chemistry-transport model GEOS-Chem. *Atmos. Environ.* **89**, 52–63 (2014).

57. Usoskin, I. G. & Kovaltsov, G. A. Production of cosmogenic ⁷Be isotope in the atmosphere: full 3-D modeling. *J. Geophys. Res. Atmos.* **113**, D12107 (2008).

58. Yeung, L. Y., Ash, J. L. & Young, E. D. Biological signatures in clumped isotopes of O₂. *Science* **348**, 431–434 (2015).

59. Koch, D. M. & Rind, D. H. Beryllium 10/beryllium 7 as a tracer of stratospheric transport. *J. Geophys. Res.* **103**, 3907–3917 (1998).

60. Field, C. V., Schmidt, G. A., Koch, D. M. & Salyk, C. Modeling production and climate-related impacts on ¹⁰Be concentration in ice cores. *J. Geophys. Res.* **111**, D15107 (2006).

61. Usoskin, I. G. et al. Short-term production and synoptic influences on atmospheric ⁷Be concentrations. *J. Geophys. Res.* **114**, D06108 (2009).

62. Usoskin, I. G., Bazilevskaya, G. A. & Kovaltsov, G. A. Solar modulation parameter for cosmic rays since 1936 reconstructed from ground-based neutron monitors and ionization chambers. *J. Geophys. Res.* **116**, A02104 (2011).

63. Clette, F. & Lefèvre, L. The new sunspot number: assembling all corrections. *Sol. Phys.* **291**, 2629–2651 (2016).

64. Fraser-Smith, A. C. Centered and eccentric geomagnetic dipoles and their poles, 1600–1985. *Rev. Geophys.* **25**, 1–16 (1987).

65. Hansen, J. et al. Climate impact of increasing atmospheric carbon dioxide. *Science* **213**, 957–966 (1981).

66. Miller, R. L. et al. CMIP5 historical simulations (1850–2012) with GISS ModelE2. *J. Adv. Model. Earth Syst.* **6**, 441–478 (2014).

67. The NCAR command language, version 6.6.2 (UCAR/NCAR/CISL/TDD, 2019); https://doi.org/10.5065/D6WD3XH5.

**Acknowledgements** This work was supported by the David and Lucile Packard Foundation Science & Engineering Fellowship and by the European Commission's Seventh Framework Programme ERC2011 under grant agreement number 291062 (ERC ICE&LASERS). We thank M. Twickler, G. Hargreaves and R. Nunn at the National Science Foundation Ice Core Facility for curating and providing ice-core samples for this study. We also thank X. Fain, A. Lemoine and G. Teste for CO₂ and CH₄ laboratory measurements at IGE; G. P. Lee and W. Sturges and his team at the University of East Anglia for halocarbon measurements on the LOCK-IN canisters; J. Freitag and C. Florian Schaller at AWI and K. Fourteau at IGE for providing a high-resolution LOCK-IN density profile; E. Le Meur for evaluating the site elevation; K. Fourteau, C. F. Schaller and J. Savarino for discussions; and A. Landais for comments on the manuscript. The LOCK-IN field and scientific programme was funded by Institut Polaire Français Paul-Emile Victor programme number 1153 and Centre National de la Recherche Scientifique INSU/LEFE programme NEVE-CLIMAT. We thank the field personnel at the LOCK-IN site: D. Colin, P. Dordhain, P. Possenti, as well as P. Godon for setting up the field logistics.

**Author contributions** L.Y.Y., L.T.M. and J.C. designed the project. J.C., P.M. and A.O. collected the firn-air samples and analysed them for trace gases and bulk isotope ratios. L.Y.Y., H.H. and A.B. performed firn and ice-core clumped-isotope measurements. L.T.M. performed additional development of the GISS-E2.1 model and conducted the global model simulations, which were then analysed by L.Y.Y. P.M. and E.W. developed the atmospheric history inversions. L.Y.Y. wrote the paper with input from all authors.

**Competing interests** The authors declare no competing interests.

**Additional information**

**Correspondence and requests for materials** should be addressed to L.Y.Y.

**Peer review information** Nature thanks Mat Evans, Jan Kaiser and the other anonymous reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at http://www.nature.com/ reprints.
Age distributions

<figure><img src="image_5.png" /></figure>

**Extended Data Fig. 1** | Age distribution estimates in firn and ice.
In green: gas age distributions in LOCK-IN firn at depths of 84.2 m (short-dashed line), 98.6 m (long-dashed line), 104.8 m depth (solid line) and 107.65 m (short-long dashed line). In red: gas age distribution

of GISP2 ice samples. In blue: gas age distributions of WAIS Divide ice samples, estimated using a diffusivity based on NEEM-EU data (solid line) and NEEM-US data (dashed line). In purple: gas age distribution of Siple Dome ice samples.
<figure><img src="image_6.png" /></figure>

**Extended Data Fig. 2 | Comparison of the most precise model (in isotopic δ notation<sup>23</sup>) and the simplified model used to include ice data<sup>45</sup>.** Left, black stars show Δ<sub>36</sub> ('D-36', in parts per thousand) data in LOCK-IN firn plotted against mean gas ages with uncertainties (±2 s.e.m., calculated using the pooled standard deviation) shown as vertical bars. Lines represent reconstructed atmospheric trends. The preferred scenario is obtained using a SCRIPPS-based O<sub>2</sub> trend (see Methods, 'Δ<sub>36</sub> for firn modelling') and is constrained by LOCK-IN firn data excluding the deepest value. The black and green solid lines show the best-guess trend obtained with the most precise and simplified

<figure><img src="image_7.png" /></figure>

models, respectively. Long-dashed lines show the uncertainty envelope. Short-dashed lines show the results of the two models when including the probably contaminated deepest data point (most precise model in green, simplified model in black). Right, δ<sup>18</sup>O<sup>18</sup>O<sub>cor</sub> data ('d18O18Ocor', in parts per thousand; see Methods, 'Δ<sub>36</sub> for firn modelling') in LOCK-IN firn against depth (symbols), compared with model results (lines). The four model simulations only differ in the very deep firn. The deepest data point at 107.65 m, which is probably contaminated (see text), is not shown in the left panel. Its mean gas age is 110 yr, corresponding to 1906 AD.
<figure><img src="image_8.png" /><figcaption>**Extended Data Fig. 3** | Effect of gravitational fractionation on isotopic ratios in the LOCK-IN firn. Shown are δ¹⁵N values of N₂ measured at LSCE (black stars); δ¹⁸O values of O₂, divided by 2, measured at LSCE</figcaption></figure>

(green crosses) and at Rice University (red circles); and δ¹⁸O¹⁸O values of O₂, divided by 4, measured at Rice University (blue circles). The black line shows the barometric slope.
<figure><img src="image_9.png" /><figcaption>Extended Data Fig. 4 | Test of gravitational fractionation corrections for δ<sup>18</sup>O<sup>18</sup>O in LOCK-IN firn. Shown are the corrections obtained using the δ<sup>15</sup>N value of N<sub>2</sub> measured at LSCE (black stars), the δ<sup>18</sup>O value of O<sub>2</sub></figcaption></figure>

measured at LSCE (green crosses), the δ<sup>18</sup>O value of O<sub>2</sub> measured at
Rice University (red circles) and the Δ<sub>36</sub> value of O<sub>2</sub> measured at Rice
University (empty blue circles).
<figure><img src="image_10.png" /></figure>

**Extended Data Fig. 5** | Comparison of firn model scenarios. Shown are results obtained with or without the deepest LOCK-IN data point, and with constant or SCRIPPS-based $O_2$ trend estimates, as well as forward firn model tests of atmospheric model scenarios. Top left, $\Delta_{36}$ data in firn and ice (LOCK-IN in green, GISP2 in red, WAIS Divide in blue, Siple Dome in purple) plotted against mean $^{18}$O$^{18}$O age, compared with atmospheric trends obtained by inverse firn/ice modelling. Shown also are the $\pm 2\sigma$-equivalent uncertainty envelope for the inverse model (long-dashed black lines) and the best-guess trends obtained using: the SCRIPPS-based $O_2$ scenario and excluding the deepest LOCK-IN data point (short-dashed black line); a constant-$O_2$ scenario and excluding the deepest LOCK-IN data point (red line); the SCRIPPS-based $O_2$ scenario and all LOCK-IN data points (dashed grey line); and a constant $O_2$ scenario and all LOCK-IN data points (blue line). Top right, $\delta^{18}O^{18}O_{cor}$ in firn and ice (defined in Methods, '$\Delta_{36}$ for firn modelling'; LOCK-IN in

<figure><img src="image_11.png" /></figure>

green, GISP2 in red, WAIS Divide in blue, Siple Dome in purple) plotted against depth, compared with model results in firn and ice using the SCRIPPS-based $O_2$ scenario. The solid lines show the simulation excluding the deepest LOCK-IN data point and the dashed lines correspond to the simulation with all data points. Bottom left, $\Delta_{36}$ data in firn and ice (same colours as in upper panels) compared with simulated profiles using the forward firn model (LOCK-IN, solid lines; WAIS Divide, long dashed lines; GISP2, short dashed lines; Siple Dome, short-long-dashed lines). Outputs shown correspond to SCRIPPS-based atmospheric concentration trends for $O_2$ and constant values for $\delta^{18}O$ and $\Delta_{36}$ (black lines); constant values for $O_2$, $\delta^{18}O$ and $\Delta_{36}$ (grey lines; results are nearly the same as the black lines); SCRIPPS-based atmospheric concentration trends for $O_2$ and constant values for $\delta^{18}O$, with the +25% box model scenario for $\Delta_{36}$ (orange lines); the +200% box model scenario for $\Delta_{36}$ (blue lines); and the +300% box model scenario for $\Delta_{36}$ (red lines).
<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 6 | Results of sensitivity tests on atmospheric trend reconstructions from the inverse firn model.** $\Delta_{36}$ data in firn and ice (stars with $\pm 2$ s.e.m. uncertainties shown as vertical bars; LOCK-IN in green, GISP2 in red, WAIS Divide in blue, Siple Dome in purple) plotted against mean $^{18}$O$^{18}$O age and compared with modelled atmospheric trends (lines). The solid black line is the preferred scenario, obtained using a SCRIPPS-based $O_2$ concentration trend and excluding the deepest LOCK-IN data point, with its uncertainty envelope shown alongside (dashed black lines). The left panel shows the simulation that includes the deepest LOCK-IN data point (red line), a simulation with the deepest LOCK-IN data point corrected (grey star) from a maximum estimate of

<figure><img src="image_13.png" /></figure>

10% surface air contamination (purple line), tests of the sensitivity to the optimal solution (grey lines; see Methods, 'Sensitivity tests on atmospheric trend reconstructions'), the simulation excluding WAIS Divide data (blue line), and simulations excluding the Siple Dome data (green solid line) or excluding the Siple Dome data and using NEEM-US-data-based diffusivity to simulate WAIS Divide firn (green dashed line). The dashed grey line shows that a straight trend with a weak slope can remain in the uncertainty envelope. The right panel shows tests of LOCK-IN firn physics parameters (green; see Methods, 'Sensitivity tests on atmospheric trend reconstructions') and tests of LOCK-IN diffusivity constrained with field data only (blue), all nearly superimposed to the preferred trend.
<figure><img src="image_14.png" /></figure>

Extended Data Fig. 7 | Comparison of Δ₃₆ values in measurements of ice-core, firn and modern air. a, Firn and ice-core Δ₃₆ values (means of replicates) plotted against mean gas age. b, Kernel-smoothed probability density distributions of bootstrap-resampled mean values of each dataset,

showing a significant (*P* < 0.002) difference between the means of the ice-core and firn (above 105 m) datasets. Uncertainties are omitted for clarity. Pooled standard deviations for each sample type are 0.03‰–0.04‰ (see Methods). C.E., common era.
<figure><img src="image_15.png" /></figure>

**Extended Data Fig. 8 | Modelled fraction of stratospheric air derived from GISS-E2.1 ⁷Be simulations.** Shown are 12-month grid-scale moving averages at the surface near each of the four polar sampling sites between 1850 AD and 2015 AD.
<figure><img src="image_16.png" /><figcaption>Extended Data Fig. 9 | Instantaneous tropospheric O₃ radiative forcing at the tropopause relative to 1850 AD. Results calculated with the GISS-E2.1 (2000s) and GEOS-Chem/MERRA2 (2005) models. Maps created using the NCAR command language⁶⁷.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

**Extended Data Fig. 10 | Solar modulation potential.** Results estimated using the sunspot number (orange line) and reconstructed from ground-based neutron monitors and ionization chambers (from Usoskin et al.⁶²; blue line).