Simulation data and software scripts used in calculus of ∆36 signature from EMAC clumped O2 isotope-inclusive model

This publication contains simulation data and software scripts for calculating quantities related to clumped oxygen isotope signature (∆36) derivation, as described in the "static" framework of (Yeung‍ et‍ al., 2016, hereinafter "Y16") and subsequently used in Yeung‍ et‍ al.‍ (2019) analysis. We provide the output of the 1950–2011 transient simulation with EMAC model (Gromov‍ et‍ al., 2019) with explicit "dynamic" simulation of ∆36 (i.e. 18O18O isotopologues undergoing transport, mixing and O(3P)-mediated isotope equilibration) to demonstrate the importance of several assumptions/simplifications involved in the Y16 calculus.

Background

Y16 framework employs stratosphere-troposphere mixing model which operates with the following variables regarded as tropospheric averages (see Eq.(2) in Y16):

Etrop       : O2 isotope exchange (equilibration) rate mediated by O(3P)

Tequil      : temperature at whish exchange (equilibration) occurs

∆36,Tequil : tropospherix mixing end-member (corresponds Tequil, i.e. equilibration in the troposphere)

∆36,strat   : end-member composition of stratospheric O2 entering the troposphere

∆36,trop    : resulting tropospheric isotope composition (under assumption that tropospheric O2 is well  mixed on annual timescales)

Using the EMAC data, we check several assumptions/simplifications involved in the Y16 calculus, which appear to underrepresent important processes in the atmosphere. These are, viz.:

- calculation of ∆36,Tequil value is calculated by merely projecting the whole-troposphere Tequil value onto Wang‍ et‍ al.‍ (2004) ∆36–Tequil curve (i.e. disregarding atmospheric mixing effects),

- use of annual instead of monthly averages of temperature and species concentrations (renders effectively lower Tequil, as O(3P) is mostly present in warmer seasons),

- use of additional "residence time" τbox (density weighting) for tropospheric integrals.

Performing calculations

Two scripts are used to calculate/retrieve "static" and "dynamic" data, which require Ferret NOAA PMEL (http://ferretop.pmel.noaa.gov/Ferret/, v7.4 or later) software.

In order to perform all operations yourself, please download all files to the same directory and run the following scripts using Ferret, e.g.:

[user@pc]/home> cd /home/user/path_where_files_are
[user@pc]/home/user/path_where_files_are> ferret
NOAA/PMEL TMAP
FERRET v7.4 (optimized)
Linux 2.6.32-696.18.7.el6.x86_64 64-bit - 04/17/18

yes? go dynamic_calc-EMAC.jnl [start_year] [end_year]
yes? go static_calc-Y16.jnl [diag_year] [start_year] [end_year]

with optional parameters:

                [start_year] and [end_year]: the span for which time series calculated

                [diag_year]: output diagnostic plots of major terms for a given year

The scripts produce several files, of which the following ones are used in subsequent analysis:

       dynamic_calc-EMAC__ATom_TR-pr4b__1950-2011_am.dat

       static_calc-Y16__ATom_TR-pr4b__1950-2011_am.nc

The naming of indicates which assumptions were kept/discarded in the below calculation of a variable , viz.:

_Y16 – implies using Y16 calculus "as is",

_MW – using only mass-weighting for tropospheric integrals (without additional density weighting) and molar (as opposed to volumetric) exchange rates,

_*_aa – using annual instead of monthly averages of temperature, O(3P) and O2 concentrations.

The trend analysis of the calculated time series is done in spreadsheet software (see static_calc-Y16_am-trends.xlsx).

Extended figures and analysis are made using Veusz (https://veusz.github.io, v3.0.1.1 or later) software (see static_calc-Y16_am.vsz), the results are presented here as well (see static_calc-Y16_am.pdf).

Results / discussion

Note the diagnostic plots highlighting differences in the calculated tropospheric integrals. One of the largest effects on the result has the choice of additional “residence time” weighting in Y16, which is assumed proportional to air density (see relative contribution of different areas in troposphere to the final averages in figure below, also in diag-rel_cont_to_trop_ave-2005.gif plot).

Our analysis of the obtained time series highlight the following critical points on the Y16 calculus:

Use of annual average temperatures and O(3P) concentrations leads to effectively lower exchange temperatures – because higher O(3P) concentrations usually coincide with higher temperatures in warmer seasons. The amount of the effect is 1.5−2 K expressed in Tequil.

In Eqs. (3) and (4): Use of Vtrop outside the integral is inappropriate, as the rate is assumed to be volumetric ([molec/cm3/s], it is followed by dV under the integral). This disregards volumetric weighting of the rate, in favour of other weightings under the integral.

Use of “residence time weighting” (τbox) is not justified – the model in Eq. (2) assumes tropospheric O2 is well-mixed (Y16). Any additional (to mass) weighting automatically implies the troposphere is not well-mixed.

Irrespective of 3), use of air density as the proxy for τbox is not justified – simply regard two boxes with very different air densities and temperatures of 0°C and 10°C, respectively. It happens that the same amount of O2 molecules underwent isotope exchange in each gridbox. What is the average equilibration temperature of molecules over both boxes? It should be 5°C, however density weighting will inappropriately shift this value towards that of the denser box.

In Eqs. (4) and (5): kexch is used in volumetric units (i.e. [molec/cm3/s], as it is multiplied by [O(3P)][O2] to obtain the exchange rate) but inappropriately weighted by air mass and density (via τbox). This is another inconsistency in weighting.

Eq. (2) is used to derive ∆36,strat that fits observed Δ36 and calculated Etrop and ∆36,Tequil. Why the same (as for troposphere) calculus is not used for the stratospheric end-member signature? Note, Eq. (2) implies that stratosphere is a well-mixed reservoir as well (and it is so on annual timescales).

Y16 calculus disregards atmospheric mixing effects. Tropospheric average of gridcell temperatures and ∆36,Tequil (calculated per cell) will not be located on Wang curve – the average ∆36,Tequil value will always be higher (about (0.065–0.067)‰ for present-day conditions) than that expected from the average Tequil.

To recap, inconsistencies in 1), 2), 3) and 5) lead to erroneous estimates of Etrop, Tequil, ∆36,Tequil and therefore implicitly ∆36,strat, which was derived as the value fitting observed ∆36 in Eq. (2) or Y16.

The major consequence is that temperatures and rates in lower-tropospheric gridcells have very high weight in the resulting averages. The correct weighting results in much larger sensitivity of ∆36,trop to upper-tropospheric O3, as shown in figure below for 2005.

Typical underestimated and correct values are listed below (based on EMAC output):

Parameter Y16 (erroneous) MW (correct) static/with mixing Difference (MW−Y16)
------------------------------------------------------------------------------------
Etrop 1.05–1.30 1.45–1.7 ~30%
[10^19 mol/a]

Tequil −(15–13) −(38–35) ~3K
[°C]

∆36,Tequil 1.99–1.96 2.34–2.29 / 2.41–2.36 ~0.06‰ / ~0.12
[‰]

∆36,strat 2.3‰ 2.8‰ ~0.5‰

T(∆36,strat) −35.7
T(O3=150 ppbv) −67.5 ~30K
[°C]

(Ranges denote change throughout 1950–2011 period)

To note, ∆36,Tequil values estimated in Y16 are (1.93–1.90)‰, which suggests even higher tropospheric equilibration temperatures (–(10.6−8.4)°C, respectively) than listed above. ∆36,strat value estimated in Y16 corresponds temperature of −35.7°C, which is ~30K higher than typical temperatures at O3=150 nmol/mol height (EMAC output). Furthermore, this value does not fit the presented observational data, when substituted into Eq. (2) of Y16.

Regarding the analysis Yeung‍ et‍ al.‍ (2019) based on Y16 calculus, we outline further critical points:

Absolute ∆36,trop value estimated using erroneous terms from Y16 calculus and inconsistent ∆36,strat is likely wrong.

Use of constant ∆36,strat (estimated for present-day conditions) implies stratospheric O3 does not change, which is not justified.

It can be shown with Eq. (2) of Y16 that different ∆36,strat values lead to different temporal trends in ∆36,trop. Hence derived ∆36,trop temporal changes are not correct as well.

If mixing is regarded, higher ∆36,Tequil will result in less strong temporal trend in estimated ∆36,trop (about ~1/2 of the static calculus value for 1950–2011, based on EMAC data).

References:

Gromov,‍ S., Röckmann,‍ T., Laskar,‍ A.‍ H.,‍ and Peethambaran,‍ R.: Modelling the abundance of 18O18O in the atmosphere and its sensitivity to temperature and O3 photochemistry, Geophys. Res. Abstr., 21, EGU2019‑16938, doi:‍ 10.5281/zenodo.3545956, 2019.

Wang,‍ Z., Schauble,‍ E.‍ A.,‍ and Eiler,‍ J.‍ M.: Equilibrium thermodynamics of multiply substituted isotopologues of molecular gases, Geochim. Cosmochim. Acta, 68, 4779−4797, doi:‍ 10.1016/j.gca.2004.05.039, 2004.

Yeung,‍ L.‍ Y., Murray,‍ L.‍ T., Ash,‍ J.‍ L., Young,‍ E.‍ D., Boering,‍ K.‍ A., Atlas,‍ E.‍ L., Schauffler,‍ S.‍ M., Lueb,‍ R.‍ A., Langenfelds,‍ R.‍ L., Krummel,‍ P.‍ B., Steele,‍ L.‍ P.,‍ and Eastham,‍ S.‍ D.: Isotopic ordering in atmospheric O2 as a tracer of ozone photochemistry and the tropical atmosphere,‍ J. Geophys. Res. Atm., 121, 12,541−512,559, doi:‍ 10.1002/2016JD025455, 2016.

Yeung,‍ L.‍ Y., Murray,‍ L.‍ T., Martinerie,‍ P., Witrant,‍ E., Hu,‍ H., Banerjee,‍ A., Orsi,‍ A.,‍ and Chappellaz,‍ J.: Isotopic constraint on the twentieth‑century increase in tropospheric ozone, Nature, 570, 224−227, doi:‍ 10.1038/s41586‑019‑1277‑1, 2019.