Trends in global tropospheric hydroxyl radical and methane lifetime since 1850 from AerChemMIP

We analyse historical (1850–2014) atmospheric hydroxyl (OH) and methane lifetime data from Coupled Model Intercomparison Project Phase 6 (CMIP6)/Aerosols and Chemistry Model Intercomparison Project (AerChemMIP) simulations. Tropospheric OH changed little from 1850 up to around 1980, then increased by around 9 % up to 2014, with an associated reduction in methane lifetime. The model-derived OH trends from 1980 to 2005 are broadly consistent with trends estimated by several studies that infer OH from inversions of methyl chloroform and associated measurements; most inversion studies indicate decreases in OH since 2005. However, the model results fall within observational uncertainty ranges. The upward trend in modelled OH since 1980 was mainly driven by changes in anthropogenic near-term climate forcer emissions (increases in anthropogenic nitrogen oxides and decreases in CO). Increases in halocarbon emissions since 1950 have made a small contribution to the increase in OH, whilst increases in aerosol-related emissions have slightly reduced OH. Halocarbon emissions have dramatically reduced the stratospheric methane lifetime by about 15 %–40 %; most previous studies assumed a fixed stratospheric lifetime. Whilst the main driver of atmospheric methane increases since 1850 is emissions of methane itself, increased ozone precursor emissions have significantly modulated (in general reduced) methane trends. Halocarbon and aerosol emissions are found to have relatively small contributions to methane trends. These experiments do not isolate the effects of climate change on OH and methane evolution; however, we calculate residual terms that are due to the combined effects of climate change and non-linear interactions between drivers. These residual terms indicate that non-linear interactions are important and differ between the two methodologies we use for quantifying OH and methane drivers. All these factors need to be considered in order to fully explain OH and methane trends since 1850; these factors will also be important for future trends.

1 Introduction

The hydroxyl radical (OH) is a highly reactive, and consequently very short-lived, component of the Earth's atmosphere that lies at the heart of atmospheric chemistry. It is often referred to as the cleansing agent of the atmosphere, as it is the main oxidant of many important trace gases, including methane (CH₄), carbon monoxide (CO), and non-methane volatile organic compounds (NMVOCs). Hydroxyl controls the removal rates of these species and hence their atmospheric residence times (e.g. Holmes et al., 2013; Turner et al., 2019). Because of this key role in determining the trace gas composition of the atmosphere, it is important to understand what controls OH's global distribution, its temporal evolution, and drivers of changes (e.g. Lawrence et al., 2001; Murray et al., 2014; Nicely et al., 2020).

The primary source of OH is from the reaction of excited oxygen atoms (O(¹D)) with water vapour; the excited oxygen originates from the photolysis of ozone (O₃) by ultraviolet (UV; wavelength <330 nm) radiation:

$$\begin{array}{}\text{(R1)}& {\displaystyle }{\mathrm{O}}_{\mathrm{3}}+h\mathit{\nu }\to {\mathrm{O}}_{\mathrm{2}}+\mathrm{O}{(}^{\mathrm{1}}\mathrm{D}),\text{(R2)}& {\displaystyle }\mathrm{O}{(}^{\mathrm{1}}\mathrm{D})+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\to \mathrm{2}\mathrm{OH}.\end{array}$$

There is rapid cycling between OH and the hydroperoxyl radical (HO₂) as well as other peroxy radicals (RO₂, e.g. CH₃O₂). For example, oxidation of CO and CH₄ (and other NMVOCs) consumes OH and generates HO₂ and RO₂:

$$\begin{array}{}\text{(R3)}& {\displaystyle }\mathrm{CO}+\mathrm{OH}\left(+{\mathrm{O}}_{\mathrm{2}}\right)\to {\mathrm{CO}}_{\mathrm{2}}+{\mathrm{HO}}_{\mathrm{2}},\text{(R4)}& {\displaystyle }{\mathrm{CH}}_{\mathrm{4}}+\mathrm{OH}+{\mathrm{O}}_{\mathrm{2}}+\mathrm{M}\to {\mathrm{CH}}_{\mathrm{3}}{\mathrm{O}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}+\mathrm{M}.\end{array}$$

Nitrogen oxides (NO and NO₂, collectively NO_x) tend to push the OH∕HO₂ ratio in the other direction through the reaction

$$\begin{array}{}\text{(R5)}& \mathrm{NO}+{\mathrm{HO}}_{\mathrm{2}}\to {\mathrm{NO}}_{\mathrm{2}}+\mathrm{OH}.\end{array}$$

However, in strongly polluted air, NO₂ becomes a dominant sink for OH through formation of nitric acid (HNO₃). Comprehensive descriptions of hydroxyl radical chemistry are given by e.g. Derwent (1996), Stone et al. (2012), and Lelieveld et al. (2016).

Levels of OH are thus influenced by ambient levels of these other species – in particular, more CH₄, CO, and NMVOCs will reduce OH, whilst more NO_x and H₂O will increase OH through ozone chemical production and the subsequent reaction of O(¹D) with H₂O (Reaction R2) to produce OH. Water vapour is a key link between physical climate and OH, but there are many others (Isaksen et al., 2009). For example, many emissions (including biogenic and anthropogenic VOCs and lightning NO_x) and chemical reactions (e.g. Reaction R4) depend on temperature and other climate variables. Photolysis rates affect OH (e.g. Reaction R1) – hence, changes in clouds and stratospheric ozone also influence OH.

The global distribution and budget of OH have been estimated by models (e.g. Spivakovsky et al., 1990; Lelieveld et al., 2016). As part of the Fifth Coupled Model Intercomparison Project (CMIP5), the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP) analysed past and future trends in simulated OH (Naik et al., 2013; Voulgarakis et al., 2013) and attributed past changes in methane to changes in anthropogenic emissions of NO_x, CH₄, CO, and NMVOCs (Stevenson et al., 2013). However, the relative influences of different processes in driving changes in global OH remain incompletely understood (e.g. Wild et al., 2020).

Evaluation of model-simulated OH requires knowledge of real-world OH. In particular, the global distribution of OH is needed to investigate quantities such as methane lifetime. Direct measurement of OH is difficult (Stone et al., 2012), and estimates of global OH can only be inferred indirectly using measurements of species such as methyl chloroform as inputs to inverse models or through assimilation of measurements of species that constrain OH, such as CH₄, CO, and NO₂, into global atmospheric chemistry models. These methods allow trends in global OH over the last few decades to be estimated, with uncertainties (see Sect. 2.2).

This study presents results from multiple transient 1850–2014 simulations performed for CMIP6 (Eyring et al., 2016) and the associated AerChemMIP (Collins et al., 2017), and it is organized as follows. Section 2 describes how CMIP6 models simulated OH, and methods used in past studies for inferring OH trends from measurements. In Sect. 3 we present pre-industrial (PI; here taken as the 1850s) and present-day (PD) zonal mean fields of modelled OH and related species, together with historical time series of global tropospheric OH and corresponding CH₄ loss rates and lifetimes, including from sensitivity experiments that isolate the effects of specific drivers. Section 4 discusses the results, comparing trends in OH from measurements and models, estimates the roles of specific drivers in the historical evolution of OH, and draws conclusions.

2 Methods

2.1 AerChemMIP CMIP6 experiments and models

We used coupled historical transient (1850–2014) model simulations from CMIP6 (Eyring et al., 2016) and various atmosphere-only historical model simulations from the associated AerChemMIP (Collins et al., 2017). Results from three global state-of-the-art Earth system models that include detailed tropospheric and stratospheric chemistry were analysed: Geophysical Fluid Dynamics Laboratory Earth System Model version 4 (GFDL-ESM4), Community Earth System Model version 2 Whole Atmosphere Community Climate Model (CESM2-WACCM), and the United Kingdom Earth System Model version 1.0, with low-resolution (N96) atmosphere and low-resolution (1^∘) ocean (UKESM1-0-LL) (Table 1).

Table 1 Basic details of the AerChemMIP models analysed in this study. For more details, see the model references.

Figure 1 Time evolution (1850–2014) of global total emissions for (a) NO_x (Tg(NO₂) yr⁻¹), (b) CO (Tg(CO) yr⁻¹), and (c) NMVOC (Tg(VOC) yr⁻¹). Orange is for biomass burning, and beige is for anthropogenic emissions (Hoesly et al., 2018). The coloured lines in the NO_x and CO panels (red for CESM2-WACCM, blue for GFDL-ESM4, and green for UKESM1-0-LL) are the total emissions for each model, including natural sources. For historical biogenic VOC emissions, see Griffiths et al. (2020, Fig. 1).

Two base historical transient experiments have been analysed: “historical” and “histSST” (Table 2). The historical runs included a fully coupled ocean and multiple ensemble members. The histSST simulations were single-member atmosphere-only runs, with monthly mean time-evolving sea surface temperatures (SSTs) and sea ice prescribed from one ensemble member of the historical simulations. Identical historical anthropogenic forcings were applied in all base runs by using prescribed long-lived greenhouse gas and halocarbon mole fractions (Meinshausen et al., 2017) as well as anthropogenic and biomass burning emissions of near-term climate forcers (NTCFs; i.e. aerosols, aerosol precursors, and ozone precursors) (van Marle et al., 2017; Hoesly et al., 2018). Emissions of NO_x, CO, and NMVOC from 1850 to 2014 are shown in Fig. 1. Natural emissions of these species were either prescribed (e.g. soil NO_x emissions, oceanic CO emissions) or internally calculated (e.g. biogenic isoprene, lightning NO_x) by embedded process-based climate-dependent schemes that differ between models (e.g. Griffiths et al., 2020; Turnock et al., 2020). Methane mole fractions were prescribed at the surface based on observations and ice core data (Meinshausen et al., 2017); away from the surface, methane was simulated by the model. However, by prescribing surface mole fractions, methane throughout the model domain was effectively prescribed.

Table 2 Number of ensemble members analysed from CMIP6 experiments in this study. All were transient 1850–2014 simulations with evolving land surfaces, trace species emissions, and greenhouse gas (GHG) mole fractions. Baseline historical runs had freely evolving oceans, whilst histSST runs were atmosphere-only with prescribed (observed) SSTs and sea ice. Sensitivity runs are based on histSST. The “-piNTCF” simulation held emissions of all NTCFs (aerosols, their precursors, and tropospheric ozone precursors) at PI levels. The “-1950HC” simulation held halocarbon mole fractions at 1950 levels (essentially PI levels). The “-piCH4” simulation held methane mole fractions at PI levels. The “-piO3” simulation held anthropogenic tropospheric ozone precursor emissions at PI levels.

^∗ UKESM1-0-LL sensitivity run for piNTCF is based on the historical (not histSST) run.

We also analysed several variants of the histSST base case with either methane mole fractions or emissions of NTCFs fixed at pre-industrial levels or halocarbon mole fractions fixed at 1950 levels. These variants allow us to estimate the roles of different drivers in changing OH (Table 2).

For some model variables we separated fields at the tropopause to provide a methane lifetime with respect to loss processes in the troposphere and stratosphere as separate values. We used World Meteorological Organization (WMO) defined tropopause pressures (the AerChemMIP variable ptp) from the models to diagnose this masking. The exact definition used is not critical, as most oxidation occurs well away from the tropopause in the tropical lower atmosphere (compare to tropospheric ozone, for which the tropopause definition is much more important; Griffiths et al., 2020).

Models diagnosed methane loss rates due to chemical destruction in each grid box – these are dominated by reaction with OH (Reaction R4) but also include other reactions, such as the reaction of methane with Cl in the stratosphere. We have used these loss rates to calculate grid box methane lifetimes. Whole-atmosphere chemical lifetimes were calculated by dividing the total methane burden by the total loss flux over the whole model domain, or just the troposphere or stratosphere, for tropospheric and stratospheric lifetimes.

We used the histSST-piNTCF simulations to diagnose the methane–OH feedback factor (Prather, 1996). These simulations held NTCF emissions at PI levels, but methane mole fractions evolved following their historical trajectory; from 1950 onwards, halocarbon mole fractions also increased. The methane–OH feedback is normally diagnosed from dedicated experiments that perturb only the methane mole fraction (Holmes, 2018), but such experiments are only available for PI conditions within AerChemMIP (e.g. Thornhill et al., 2020a). The methane–OH feedback factor, f, was calculated as follows:

$$\begin{array}{}\text{(1)}& f={\left(\mathrm{1}-\left({\displaystyle \frac{\ln \left(\frac{{\mathit{\tau }}_{\mathrm{1930}-\mathrm{1960}}}{{\mathit{\tau }}_{\mathrm{1850}}}\right)}{\ln \left(\frac{{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{\mathrm{1930}-\mathrm{1960}}}{{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{\mathrm{1850}}}\right)}}\right)\right)}^{-\mathrm{1}},\end{array}$$

where τ is the total methane lifetime, additionally including a soil sink; CH₄ is taken to have a lifetime with respect to soil uptake of 160 years based on results for the 2000s from Spahni et al. (2011), Ito et al. (2012), Kirschke et al. (2013), and Tian et al. (2015), as summarized in Tian et al. (2016); the tropospheric Cl sink for methane is neglected here (Allan et al., 2007; Hossaini et al., 2016; Sherwen et al., 2016; Wang et al., 2019; Gromov et al., 2018; Strode et al., 2020); and [CH₄] is the global mean methane mole fraction for a particular year (or range of years) of the histSST-piNTCF simulation. The reference year is 1850, the first year of the simulation. We took average values between 1930 and 1960 to give the most reliable estimate of f, as this is after a sufficiently large methane perturbation has built up but before halocarbons interfere with the results in these simulations (see Sect. 3.3).

We used each model's feedback factor to calculate equilibrium PD methane mole fractions ([CH₄]_eq) for each sensitivity run using the diagnosed total methane lifetimes from these experiments. The equilibrium methane mole fraction is the methane mole fraction that would have been reached if methane mole fractions had not been prescribed in these runs, but rather that methane emissions had been applied, allowing methane mole fractions to evolve freely (e.g. Fiore et al., 2009; Stevenson et al., 2013):

$$\begin{array}{}\text{(2)}& [{\mathrm{CH}}_{\mathrm{4}}{]}_{\text{eq}}=[{\mathrm{CH}}_{\mathrm{4}}{]}_{\text{ref}}{\left({\displaystyle \frac{{\mathit{\tau }}_{\text{PD}}}{{\mathit{\tau }}_{\text{ref}}}}\right)}^f,\end{array}$$

where [CH₄]_ref is the prescribed methane mole fraction in the run, and τ_ref is the total methane lifetime in the histSST base experiment either for PD or, in the case of histSST-piCH4, for PI. We illustrate this with two examples: (i) in the histSST-piNTCF case, the equilibrium value is the PD methane mole fraction that would have been reached if all NTCF emissions been held at PI levels, whilst CH₄ emissions had followed their historical evolution; and (ii) in the histSST-piCH4 case, the equilibrium value is the PD methane mole fraction that would have been reached if CH₄ emissions had been held at PI levels, but all other emissions followed their historical evolution. This allows us to clarify modelled influences on CH₄ from specific emissions.

2.2 Inferred OH from measurements

Tropospheric OH has a chemical lifetime of less than a second or so, reflecting its high reactivity and making direct measurement difficult and impractical for constraining global OH distributions (e.g. Stone et al., 2012). Instead, tropospheric mean OH and its variability have traditionally been inferred from measurements of trace gases with lifetimes longer than the timescale of tropospheric mixing and whose primary loss is via reaction with OH. If emissions are well-known, then observed changes in atmospheric abundance may be related, via inverse methods, to variations in OH. To date, the favoured proxy for estimating OH has been from measurements of methyl chloroform (1,1,1-trichloroethane; CH₃CCl₃; MCF), a synthetic industrial solvent that was banned in the late 1980s as a stratospheric-ozone-depleting substance (Lovelock, 1977; Singh, 1977; Spivakovsky et al., 1990, 2000; Montzka et al., 2000; Prinn et al., 2001). The inversions have typically spatially represented the global atmosphere as a few boxes.

The earliest MCF inversions predicted relatively large OH variability, reflecting high sensitivity to the uncertainty in residual MCF emissions (Bousquet et al., 2005; Prinn et al., 2005, 2001; Krol and Lelieveld, 2003; Krol et al., 2003). However, Montzka et al. (2011) demonstrated that by the late 1990s, residual emissions had declined sufficiently so as to be a minor source of uncertainty and that OH varied by at most a few percent in year-to-year variability. More recently, multi-box models have been used with Bayesian inverse methods to simultaneously optimize OH and MCF emissions to match MCF observations from NOAA and the Advanced Global Atmospheric Gases Experiment (AGAGE) networks, as well as multi-species inversions including methane and methane isotopologues as additional constraints (Rigby et al., 2017; Turner et al., 2017). Naus et al. (2019) further investigated the inversion methods used by Rigby et al. (2017) and Turner et al. (2017), confirming that the derivation of OH from MCF and CH₄ is quite poorly constrained, and found that OH trends with a range of different magnitudes and signs were consistent with the available data.

Some inversion studies have used models with greater spatial resolution. McNorton et al. (2016) performed inverse modelling using a 3-D chemistry-transport model (CTM) constrained by MCF data and found that OH increases contributed significantly to the slowdown in the global CH₄ growth rate between 1999 and 2006 and that the post-2007 increases in CH₄ growth rate were poorly simulated if OH variations were ignored. McNorton et al. (2018) extended this work with further constraints from Greenhouse Gases Observing Satellite (GOSAT) CH₄ and δ¹³CH₄ and found that the post-2007 CH₄ growth rate surge was most likely due to a combination of a decrease (−1.8 ± 0.4 %) in global OH and increases in CH₄ emissions.

Figure 2 (a) Time evolution of global annual mean tropospheric OH (1850–2014) expressed as a percentage anomaly relative to the 1998–2007 mean (and ensemble spreads where available) for GFDL-ESM4 (blue), UKESM1-0-LL (green), CESM2-WACCM (red), and the multi-model mean (black). (b) Observation-based inversions of global annual mean tropospheric OH for 1980–2015 from Montzka et al. (2011), Rigby et al. (2017), Turner et al. (2017), and Nicely et al. (2018), including ± 1 standard deviation uncertainties for the results from Rigby et al. (2017), with model results from (a) overlain.

These inversion studies generally find OH to have increased from the late 1980s until the mid-2000s when OH then began to decline (Fig. 2b). However, most inversion studies also found that solutions exist within the uncertainty of the system when OH was held constant and only emissions of the reactants were allowed to be optimized. In contrast, Nicely et al. (2018) empirically derived a historic global mean OH reconstruction by taking a baseline forward OH simulation from the National Aeronautics and Space Administration (NASA) Global Modeling Initiative (GMI) CTM driven by assimilated meteorology since 1980 and adjusting it based on box-model-derived relationships of OH responses to changes in observable parameters such as total ozone columns from satellites. The empirically derived OH reconstruction (also shown in Fig. 2b) was found to be relatively invariant when compared to other MCF inversions over the past few decades, which the study suggested reflected chemical buffering of the many competing factors that can influence OH.

Several studies have investigated the constraints imposed on OH by species other than MCF and CH₄. Gaubert et al. (2017) assimilated time series of global-scale satellite CO measurements from the Measurements of Pollution In The Troposphere (MOPITT) project into a global model and found a decrease in the global CO burden of ∼ 20 % over the period 2002–2013. Associated with this decrease in CO was an 8 % shortening of the methane lifetime and a 7 % increase in OH. Nguyen et al. (2020) also found that decreasing global CO concentrations since the 2000s have important influences on CH₄ flux inversion results because of the strong chemical coupling between CO, CH₄, and OH. Assimilation of satellite CO, NO_x, and O₃ data (e.g. Miyazaki et al., 2015, 2017; Miyazaki and Bowman, 2017; Gaubert et al., 2017) demonstrates that OH is sensitive to all these species and that data assimilation improves simulation of the hemispheric ratio of OH (Patra et al., 2014).

Collectively, these earlier studies have shown that OH is influenced by CO, NO₂, O₃, and CH₄. To date, studies have used subsets of the available observational data (i.e. one or more of MCF, CH₄, δ¹³CH₄, CO, NO₂, and O₃) to constrain OH but not yet all available relevant data. The OH trends derived from several of these studies, including the uncertainty estimates from Rigby et al. (2017), are summarized in Fig. 2b.

3 Results

3.1 Pre-industrial to present-day base simulations

Figure 2a shows time series (1850–2014) of global annual mean tropospheric OH burden expressed as a percentage anomaly relative to the 1998–2007 mean value for the three models. This shows typical inter-annual variability in global OH of about ±2 %–3 %, a small decrease (about −3 %) in OH from 1850 up to 1910, then a similar magnitude increase up to the 1980s. From the 1980s to 2014, the models show a strong increase in OH of about 9 %. All three models show comparable behaviour. We find very similar results between the fully coupled (historical) and the atmosphere-only (histSST) experiments (not shown). This confirms that it is valid to directly compare and analyse together the results from these two experimental set-ups.

Figure 2b shows several estimates of global tropospheric OH trends over the period 1980–2014 inferred from observations (as described in Sect. 2.2), including an uncertainty range from Rigby et al. (2017). The inferred trends from different inversion methods show quite a wide range but are generally upwards from 1980 to 2005, in broad agreement with the AerChemMIP models. However, from 2005 onwards, the inversions generally indicate downwards trends, whereas the models suggest a continued slight upwards trend. The 1980–2015 model global OH trends are almost always within the ±1 standard deviation uncertainty range from Rigby et al. (2017), although they are close to the lower end of the range in 1980 and just beyond its upper end in 2015.

Figure 3 First column: zonal mean (latitude–pressure, hPa) cross sections for OH concentration (10⁵ molec. cm⁻³) averaged over the period 2005–2014 (PD) for the histSST simulations. Rows show results for CESM2-WACCM, GFDL-ESM4, and UKESM1-0-LL. Solid lines indicate the tropopause (PD in black; other in green). Other panels show differences (%) between experiments. Second column: histSST PD minus PI (1850–1859 mean). Third column: histSST minus histSST-piNTCF for PD. Fourth column: histSST minus histSST-1950HC for PD.

Figure 3 shows present-day (PD; 2005–2014 decadal mean) zonal mean OH concentrations for the three models. The vertical coordinate is pressure, and the zonal mean WMO tropopause is indicated. All models show high OH values between 30^∘ S and 30^∘ N in the lower to middle troposphere, with larger values in the Northern Hemisphere (NH).

Figure 3 also shows changes in OH from pre-industrial (PI; 1850–1859 decadal mean) to PD expressed as the percentage change relative to PD. This reveals local increases of over 50 % in zonal mean tropospheric OH, in particular over polluted NH mid-latitudes, but also a local decrease of over 10 % in the Southern Hemisphere (SH) middle to upper troposphere at around 20^∘ S. The PD–PI figures also show both the PD and PI tropopauses and indicate insignificant changes in tropopause height over the historical era.

Figure 4 The same layout as Fig. 3, but for CH₄ lifetime (years).

Figure 4 shows the zonal mean distribution of local methane lifetime, which ranges from about 2.5 years in the tropical lower troposphere to > 20 years in colder, drier high latitudes and in the vicinity of the tropopause. Short lifetimes also occur higher in the stratosphere but do not significantly contribute to the whole-atmosphere chemical lifetime due to the low air densities at high altitudes. The multi-model mean whole-atmosphere PD chemical lifetime in histSST is 8.4 ± 0.3 years, lower than the mean PI lifetime of 9.5 ± 0.5 years (lifetimes for individual models are given in Table 3; the ranges are the standard deviations across the models). These values compare to a whole-atmosphere methane lifetime for 2010 (mean ± 1 standard deviation) of 9.1 ± 0.9 years (Prather et al., 2012), as used by the Intergovernmental Panel on Climate Change (IPCC; Myhre et al., 2013). Lifetimes have fallen since the PI, mainly reflecting increases in OH.

Table 3 Whole-atmosphere methane chemical (not including soil sink) lifetimes (years). PI refers to the 1850–1859 mean; PD refers to the 2005–2014 mean. Uncertainties are ±1 standard deviation based on the range of annual values.

^∗ UKESM1-0-LL methane lifetime for piNTCF is based on the historical (not histSST) run.

3.2 Historical sensitivity simulations

The drivers of these changes in OH and methane lifetime were explored further using results from sensitivity experiments based on the histSST simulations. These kept anthropogenic emissions or mole fractions of particular species, or groups of species, at their PI or 1950 levels (Table 2). Figure 3 shows how zonal mean OH in the models responded to fixing NTCF emissions at PI levels and halocarbon mole fractions at 1950 levels. The panels in Fig. 3 shows percentage changes in OH relative to the PD histSST base case. Figure 4 shows percentage changes in methane lifetime.

We define the annual tropospheric OH burden anomaly in the base histSST simulations at time t (ΔOH_Base(t)) as the percentage change in OH since PI (1850–1859):

$$\begin{array}{}\text{(3a)}& {\mathrm{\Delta }\mathrm{OH}}_{\text{Base}}\left(t\right)=\mathrm{100}\mathit{\%}\times \left({\displaystyle \frac{{\text{OH}}_{\text{histSST}}\left(t\right)-{\text{OH}}_{\text{histSST}}\left(\text{PI}\right)}{{\text{OH}}_{\text{histSST}}\left(\text{PI}\right)}}\right),\end{array}$$

or, for clarity, dropping the (t) and substituting OH_PI for OH_histSST(PI):

$$\begin{array}{}\text{(3b)}& {\mathrm{\Delta }\mathrm{OH}}_{\text{Base}}=\mathrm{100}\mathit{\%}\times \left({\displaystyle \frac{{\mathrm{OH}}_{\text{histSST}}-{\mathrm{OH}}_{\text{PI}}}{{\mathrm{OH}}_{\text{PI}}}}\right).\end{array}$$

We then use each sensitivity run to isolate the contributions to this overall OH anomaly from changes in CH₄ mole fraction, NTCF emissions, halocarbon mole fraction, and O₃ precursor emissions since 1850:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathrm{\Delta }\mathrm{OH}}_{{\mathrm{CH}}_{\mathrm{4}}}=\mathrm{100}\mathit{\%}\times \left({\displaystyle \frac{{\mathrm{OH}}_{\text{histSST}}-{\mathrm{OH}}_{\text{histSST\_piCH4}}}{{\text{OH}}_{\text{PI}}}}\right),\text{(5)}& {\displaystyle }{\mathrm{\Delta }\mathrm{OH}}_{\text{NTCF}}=\mathrm{100}\mathit{\%}\times \left({\displaystyle \frac{{\mathrm{OH}}_{\text{histSST}}-{\mathrm{OH}}_{\text{histSST\_piNTCF}}}{{\mathrm{OH}}_{\text{PI}}}}\right),\text{(6)}& {\displaystyle }{\mathrm{\Delta }\mathrm{OH}}_{\text{HC}}=\mathrm{100}\mathit{\%}\times \left({\displaystyle \frac{{\mathrm{OH}}_{\text{histSST}}-{\text{OH}}_{\text{histSST\_1950HC}}}{{\text{OH}}_{\text{PI}}}}\right),\text{(7)}& {\displaystyle }{\mathrm{\Delta }\mathrm{OH}}_{\mathrm{O}\mathrm{3}}=\mathrm{100}\mathit{\%}\times \left({\displaystyle \frac{{\mathrm{OH}}_{\text{histSST}}-{\text{OH}}_{\text{histSST\_piO3}}}{{\text{OH}}_{\text{PI}}}}\right).\end{array}$$

Since the ΔOH_NTCF and ΔOH${}_{{\mathrm{O}}_{\mathrm{3}}}$ anomalies only differ in that the former includes the effects of aerosols, then assuming that the impacts of aerosols and O₃ precursors on OH do not interact with each other, we can also isolate the contribution from changes in aerosols to the overall OH anomaly:

$$\begin{array}{}\text{(8)}& {\mathrm{\Delta }\mathrm{OH}}_{\text{aerosol}}={\mathrm{\Delta }\mathrm{OH}}_{\text{NTCF}}-{\mathrm{\Delta }\mathrm{OH}}_{{\mathrm{O}}_{\mathrm{3}}}.\end{array}$$

In addition, we can calculate a “residual” contribution, i.e. the component of the overall OH anomaly that is left after linearly adding all the other components:

$$\begin{array}{}\text{(9)}& {\mathrm{\Delta }\mathrm{OH}}_{\text{residual}}={\mathrm{\Delta }\mathrm{OH}}_{\text{Base}}-{\mathrm{\Delta }\mathrm{OH}}_{{\mathrm{CH}}_{\mathrm{4}}}-{\mathrm{\Delta }\mathrm{OH}}_{\text{NTCF}}-{\mathrm{\Delta }\mathrm{OH}}_{\text{HC}}.\end{array}$$

This residual component represents the contribution of climate change to the OH anomaly, along with any contributions from non-linear interactions between components. Non-linearities may arise, for example, because the response of OH to changes in CH₄ is likely to differ depending on whether NTCFs, such as NO_x, are at PI or PD levels. Such interactions are not isolated by our methodology, and it is unclear whether the climate change signal or the effects of non-linearities dominate this residual term.

Figure 5 The histSST (base, in black) tropospheric OH anomaly (percent change relative to PI) for each year (see Eq. 3) for the three models: (a) CESM2-WACCM, (b) GFDL-ESM4, and (c) UKESM1-0-LL. The coloured lines represent the contributions to this OH anomaly due to changes since 1850 in CH₄ mole fraction, NTCF emissions, halocarbon mole fraction, O₃ precursor emissions, and aerosols (see Eqs. 4–8; only NTCF and HC experiments from CESM2-WACCM). The residual curve (see Eq. 9) is the extra contribution required after linearly adding the curves for CH₄, NTCF, and HC that is needed to reproduce the base anomaly.

Figure 5 shows time series of how the base OH anomaly (Eq. 3) evolves, together with each of the components (Eqs. 4–9) that contribute to the base anomaly. Figure 6 compares the magnitudes of these various drivers of OH changes over two time periods: 1850–1980 and 1850–2010. Figure 7 shows the evolution of whole-atmosphere methane lifetime for the base histSST runs and each sensitivity run. Figure 7 also separates the methane lifetime into its tropospheric and stratospheric components.

Figure 6 Summary of drivers of OH changes (%) relative to 1850 for the three models and their multi-model mean over (a) 1850–1980 and (b) 1850–2010. We have used decadal means: 1850 refers to 1850–1859, 1980 is 1975–1984, and 2010 is 2005–2014. The shaded areas show the split of the NTCF signal (green) into ozone precursors (blue) and aerosols (brown) when models have performed both the histSST-piNTCF and histSST-piO3 experiments. The residual values (pale blue) are the differences between the total change (black, from the histSST simulations) and the sum of the changes from CH₄ (red), NTCF, and halocarbons (purple). We interpret the residual terms as being due to climate change, in addition to any non-linear interactions between forcings.

Figures 5 and 6 show that the evolution of OH has been mainly controlled by the balance between the growth of methane, which has acted to reduce OH by over 20 %, and the changes in NTCF emissions (and in UKESM1-0-LL, the residual term), which have tended to increase OH. Because these opposing drivers have similar magnitudes, small mismatches between them are key, and the other minor drivers can also be important contributors to the overall trend in OH.

The impact of increases in NTCF emissions since 1850 up to PD was to generally increase tropospheric OH by 10 %–50 % in the zonal mean (Fig. 3) and 13 %–22 % across the whole troposphere (Figs. 5 and 6); this mainly reflects the dominant role of NO_x increases, whose impact overwhelms the impacts of increasing CO (up to ∼ 1990) and NMVOC emissions, which will have tended to reduce OH. Since about 1990, global CO emissions have been reduced (Fig. 1), also contributing to the increase in OH. The overall impact of changed emissions of NTCFs has been to reduce the methane lifetime (Figs. 4 and 7, Table 3). This is mainly driven by increases in NO_x emissions. The structure seen in the zonal mean PD–PI change in OH (Fig. 3, column 2) can be largely explained by the change in NTCF emissions (Fig. 3, column 3), with the effects of methane mole fraction increases superimposed. Note that increasing methane also increases CO, and both these reduce OH. We are unable to isolate the effects of CO in our experiments.

Emissions of halocarbons since 1950 have led to polar stratospheric ozone depletion, mainly in the SH. This has increased stratospheric OH levels but also increased tropospheric OH due to increased penetration of ultraviolet (UV) radiation and consequently higher photolysis rates (Fig. 3). The overall impact on tropospheric OH and methane lifetime is comparatively small (Figs. 4–7, Table 3), but the impact on methane lifetime in the stratosphere has been dramatic, reducing it from ∼ 170 to ∼ 140 years (CESM2-WACCM), from ∼ 140 to ∼ 80 years (GFDL-ESM4), and from ∼ 190 to 145 years (UKESM1-0-LL) (Fig. 7). These changes are mainly driven by changes in stratospheric Cl. These values can be compared to an assumed constant value for the lifetime of methane with respect to stratospheric chemical destruction of 120 (±20 %) years in the IPCC Fifth Assessment Report (Prather et al., 2012).

Figure 7 Time evolution (1850–2014) of CH₄ lifetime (years) for (a–c) CESM2-WACCM, (d–f) GFDL-ESM4, and (g–i) UKESM-0-LL averaged over the whole atmosphere (a, d, g), the troposphere (b, e, h), and the stratosphere (c, f, i). Colours refer to different model experiments, as indicated in (f). For UKESM-0-LL, we used historical-piNTCF as histSST-piNTCF was not available.

The effects of increased emissions of aerosols and aerosol precursors can be diagnosed by differencing the piO3 and piNTCF simulations (Eq. 8). Aerosols have slightly reduced OH (Figs. 5 and 6) and lengthened the methane lifetime (Fig. 7), but the effect is small in magnitude compared to most other effects.

For the two models able to diagnose the residual term, they both suggest a positive impact on OH although by variable amounts (6 %–13 %), with a larger residual term in UKESM1-0-LL. We suggest this term may reflect increases in humidity associated with climate change and an increase in the primary OH production flux (Reaction R2). However, exactly what the residual terms represent remains uncertain.

3.3 Contribution of OH drivers to PI–PD changes in methane

Figure 7 shows values for the methane–OH feedback factor (from a modified version Eq. (1) using values for individual years rather than 1930–1960) calculated for every year in the histSST-piNTCF simulations. In the first few decades, the methane changes are small and the variability of the methane lifetime yields large fluctuations in f. Beyond about 1960, changes in halocarbon mole fractions mean that the values of f are unreliable. We therefore use the average value over the time period 1930–1960 as our best estimate of the feedback factor. This yields a value of 1.25 for CESM2-WACCM and 1.23 for GFDL-ESM4. Thornhill et al. (2020a) find values of f from the piClim simulations of 1.30 for GFDL-ESM4 and 1.32 for UKESM1-0-LL. The values derived using Eq. (1) are probably slightly smaller because the histSST-piNTCF runs also include increases in temperature and humidity. These values are similar to the range of values found in previous studies: 1.23–1.35 (Stevenson et al., 2013; six models), 1.19–1.28 (Voulgarakis et al., 2013; two models, year 2000 conditions), and 1.33–1.45 (Prather et al., 2001; seven models). Using the values of f for 1930–1960 (Fig. 7) for CESM2-WACCM and GFDL-ESM4, the value of 1.32 for UKESM1-0-LL (Thornhill et al., 2020a), and the lifetimes presented in Table 3, we calculate equilibrium PD methane mole fractions for all sensitivity experiments (Table 4).

Table 4 Equilibrium PD global mean methane mole fractions (ppbv) inferred from PD methane lifetimes from the sensitivity experiments. Also shown are percentage changes compared to the observed PD value (1794 ppbv) or, for the piCH4 case, the observed PI value (808 ppbv).

Observed PI and PD methane levels are 808 and 1794 ppb, respectively. Holding NTCFs at PI levels increases PD methane by 16 %–33 %. This is more intuitively interpreted in terms of the impact of the increased emissions of NTCFs: they have tended to reduce PD methane by this amount. Similarly, the impact of halocarbon emissions has been to reduce PD methane by 7 %–15%.

Taking the average of results from the GFDL-ESM4 and UKESM1-0-LL models (that have values for all categories), holding methane emissions at PI levels would have led to PD methane levels of 516 ppbv, which is 36 % (292 ppbv) lower than PI mole fractions. Hence, the net impact of increasing methane emissions has been to increase methane mole fractions from 516 to 1794 ppbv, an increase of 1278 ppbv $\left(\mathrm{\Delta }{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{{\mathrm{eq}}_{{\mathrm{CH}}_{\mathrm{4}}}}\right)$. This increase is 30 % larger than the simple observed PI to PD increase in methane (986 ppbv, $\mathrm{\Delta }{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{{\mathrm{eq}}_{\mathrm{obs}}}$). The impact of NTCF emissions was to reduce PD methane by 480 ppbv ($\mathrm{\Delta }{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{{\mathrm{eq}}_{\mathrm{NTCF}}}$), whilst increases in halocarbon emissions reduced PD methane by 149 ppbv ($\mathrm{\Delta }{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{{\mathrm{eq}}_{\mathrm{HC}}}$). These diagnosed contributions do not linearly add up to give the observed total; there is a residual term, as also found when attributing the OH changes to drivers. Following a similar format to Eq. (9), we can diagnose this residual term:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& \mathrm{\Delta }{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{{\mathrm{eq}}_{\mathrm{residual}}}=\mathrm{\Delta }{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{{\mathrm{eq}}_{\mathrm{obs}}}-\mathrm{\Delta }{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{{\mathrm{eq}}_{\mathrm{CH}\mathrm{4}}}\\ & -\mathrm{\Delta }{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{{\mathrm{eq}}_{\mathrm{NTCF}}}-\mathrm{\Delta }{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{{\mathrm{eq}}_{\mathrm{HC}}}\\ & =\mathrm{986}-\mathrm{1278}+\mathrm{480}+\mathrm{149}=\mathrm{337}\text{ppbv}.\end{array}\end{array}$$

We tentatively attributed the OH residual term to climate change impacts, as the residual OH increase could physically be linked to water vapour increases. However, the residual change in equilibrium methane is positive, whilst it would be expected to be negative in order to match the positive residual OH term. The other attributions for OH and equilibrium methane are more well-behaved and consistent. This suggests that non-linear interactions between drivers are important and differ in strength between our attribution methodologies for OH and methane. This means that perfect quantitative attribution cannot be achieved, and attribution of the residual term to climate change effects is rather uncertain. Nevertheless, the magnitudes of these attribution terms are useful qualitative indicators of the relative importance of different drivers of changes in OH and methane lifetime.

4 Discussion and conclusions

Modelled OH trends presented in this study are from state-of-the-art Earth system models driven by CMIP6 historical forcings, including observed trends in CH₄ and halocarbon mole fractions. The latter drive stratospheric ozone depletion in the models, which strongly influences tropospheric UV levels and hence photolysis rates. Apart from CH₄, all other reactive species that control OH (e.g. CO, O₃, NO₂, and H₂O) freely evolve in the simulations in response to prescribed CMIP6 emissions and simulated climate. These model simulations of OH are very important for understanding past trends and projecting future trends in CH₄.

The base model simulations all show similar historical trends in global OH, with relative stability from 1850 up to 1980, followed by strong (9 %) increases up to the present day (Fig. 2). The earlier stability is in good agreement with previous studies (e.g. Naik et al., 2013). The increase from 1980 to 2005 is broadly consistent with several studies that use MCF and other species to reconstruct OH trends from observations; however, since 2005 most of these reconstructions indicate a decrease in OH, whereas our models indicate a continued increase (Fig. 2b). However, these reconstructions show a wide range of trends, and our modelled trends fall just about within the uncertainty range estimated by Rigby et al. (2017). The magnitudes of the model's recent increases are similar to results from Gaubert et al. (2017), who assimilated satellite-derived trends in CO since 2002 into an Earth system model. Several OH inversions have used multiple observational datasets (Miyazaki et al., 2015; McNorton et al., 2018), and as the time series of observations, particularly satellite data, lengthens, uncertainties on real-world OH trends will hopefully be reduced, providing stronger constraints for models.

Figure 8 Calculated values for the methane–OH feedback factor (f) from the histSST-piNTCF experiments for CESM2-WACCM and GFDL-ESM4. Mean and standard deviation values for 1930–1960 (shaded time period) are shown.

We attempted to quantify the component drivers of the changes in OH using a series of idealized model sensitivity experiments. These experiments exhibit relatively consistent OH responses across the models (Figs. 5 and 6) and show that the evolution of methane and ozone precursor emissions have strongly influenced OH trends. Halocarbon and aerosol-related emissions have had relatively small impacts. We also diagnose a residual component that represents the impact of climate change and non-linear interactions between drivers. Other studies have indicated that climate variations and change influence OH (e.g. Naik et al., 2013; Murray et al., 2014; Turner et al., 2018). The modelled increase in OH since 1980 is because the influence of NTCF emissions, together with this residual term, outweighs the effects of increasing CH₄ (Fig. 6). These experiments did not separate the effects of different ozone precursors, but these have been explored in previous studies (Stevenson et al., 2013; Holmes et al., 2013), wherein increases in anthropogenic NO_x emissions have been found to be the main NTCF driver of OH increases. Recent reductions in anthropogenic CO emissions (Fig. 1) are clearly also important (Gaubert et al., 2017), but our experiments are unable to separate the relative impacts of these two species.

The trends in OH are associated with trends in methane lifetime (Fig. 7), and we have used these to estimate the influence of individual drivers on the methane mole fraction by calculating equilibrium methane levels from the changes in lifetime (Table 4). Drivers that increase OH lead to reductions in methane lifetime and equilibrium methane. The residual component for OH is positive and may mainly physically represent the rise of water vapour associated with climate warming. This finding is broadly consistent with results presented by Thornhill et al. (2020b) of the negative impacts on methane lifetime found in 4×CO₂ experiments (see their Table 15). However, the residual component we diagnosed from changes in equilibrium methane is also positive, which suggests that non-linear interactions show different impacts in our two methodologies that diagnose residual effects and that the residual term may not be a good indicator of climate change effects alone. These results indicate that methodologies to isolate drivers of OH and methane changes need careful interpretation, as non-linearities (i.e. couplings between drivers) appear to be important.

Although halocarbon emissions have had quite small effects on the whole-atmosphere methane lifetime, they have had dramatic impacts on methane's stratospheric chemistry, and its lifetime may have been reduced by up to about 40 % between 1960 and 1990 (Fig. 7). Previous studies have generally assumed a fixed stratospheric lifetime for methane (e.g. Prather et al., 2012).

All these factors need to be included in holistic assessments of OH and methane change. The CMIP6/AerChemMIP model simulations contain many useful diagnostics that will allow us to better understand the drivers of atmospheric OH and methane trends. This study represents a very preliminary initial analysis of this rich multi-model, multi-experiment dataset.