Atmospheric composition and climate impacts of a future hydrogen economy

Hydrogen is expected to play a key role in the global energy transition to net zero emissions in many scenarios. However, fugitive emissions of hydrogen into the atmosphere during its production, storage, distribution and use could reduce the climate benefit and also have implications for air quality. Here, we explore the atmospheric composition and climate impacts of increases in atmospheric hydrogen abundance using the UK Earth System Model (UKESM1) chemistry–climate model. Increases in hydrogen result in increases in methane, tropospheric ozone and stratospheric water vapour, resulting in a positive radiative forcing. However, some of the impacts of hydrogen leakage are partially offset by potential reductions in emissions of methane, carbon monoxide, nitrogen oxides and volatile organic compounds from the consumption of fossil fuels. We derive a refined methodology for determining indirect global warming potentials (GWPs) from parameters derived from steady-state simulations, which is applicable to both shorter-lived species and those with intermediate and longer lifetimes, such as hydrogen. Using this methodology, we determine a 100-year global warming potential for hydrogen of 12 ± 6. Based on this GWP and hydrogen leakage rates of 1 % and 10 %, we find that hydrogen leakage offsets approximately 0.4 % and 4 % respectively of total equivalent CO₂ emission reductions in our global hydrogen economy scenario. To maximise the benefit of hydrogen as an energy source, emissions associated with hydrogen leakage and emissions of the ozone precursor gases need to be minimised.

1 Introduction

The adoption of low-carbon hydrogen (H₂) as an energy source could lead to substantial reductions in carbon dioxide emissions and a significant climate benefit. However, increasing levels of hydrogen in the atmosphere as a result of fugitive emissions from a hydrogen economy will affect atmospheric composition (e.g. Prather, 2003). It is important that the atmospheric implications of potential changes in hydrogen emissions are investigated in detail before the implementation of widespread hydrogen use.

The hydrogen abundance in the atmosphere increased over the twentieth century (Patterson et al., 2021) and today has a mixing ratio of around 530 parts per billion (ppb), with a small interhemispheric gradient (±20 ppb) (Novelli et al., 1999; Patterson et al., 2021). Fossil fuel combustion and biomass burning account for approximately 50 % of the current total global hydrogen source, with the remainder arising from the oxidation of methane (CH₄) and volatile organic compounds (VOCs) in the atmosphere (e.g. Ehhalt and Rohrer, 2009; Pieterse et al., 2013; Grant et al., 2010). Hydrogen removal is dominated primarily by uptake to soils and reaction with hydroxyl radicals (OH). OH is the main atmospheric oxidising agent, and if hydrogen emissions were to increase, subsequent changes in the OH concentration could alter the lifetimes of important atmospheric greenhouse gases. Therefore, whilst hydrogen itself is not radiatively active, it can act as an indirect greenhouse gas.

An increase in the tropospheric concentration of hydrogen reduces the availability of OH via Reaction (R1). Reductions in tropospheric OH will result in increases in the atmospheric lifetime of CH₄ and its abundance (Derwent et al., 2001; Schultz et al., 2003; Warwick et al., 2004; Derwent et al., 2020; Field and Derwent, 2021) since the major atmospheric sink of CH₄ is through reaction with OH (Reaction R2).

$$\begin{array}{}\text{(R1)}& {\displaystyle }{\mathrm{H}}_{\mathrm{2}}+\mathrm{OH}\to {\mathrm{H}}_{\mathrm{2}}\mathrm{O}+\mathrm{H}\text{(R2)}& {\displaystyle }{\mathrm{CH}}_{\mathrm{4}}+\mathrm{OH}\to {\mathrm{CH}}_{\mathrm{3}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\end{array}$$

The oxidation of both hydrogen and CH₄ in the troposphere can lead to the generation of ozone, which impacts both climate and air quality (e.g. Archibald et al., 2020a; Schultz et al., 2003). In the stratosphere, oxidation of hydrogen and methane will lead to increases in stratospheric water vapour with potential implications for climate and stratospheric ozone (e.g. Tromp et al., 2003; Warwick et al., 2004). Increases in atmospheric hydrogen could also influence stratospheric ozone recovery through the production of HO_x radicals (i.e. OH + HO₂), which are involved in ozone destruction cycles.

The extent to which future changes in hydrogen might affect atmospheric composition and climate will depend upon the level of hydrogen leakage and the ultimate size of a future hydrogen economy. In addition, emission reductions in species currently emitted by the production and consumption of fossil fuels, including methane, CO, NO_x (i.e. NO + NO₂) and volatile organic compounds (VOCs), will also induce feedbacks on atmospheric composition (e.g. Jacobson, 2008). There is significant uncertainty in the size of a future hydrogen economy, and the leakage rates and emission reductions in other species are uncertain, depending on the forms of new technology implemented (e.g. Frazer-Nash Consultancy, 2022; Lewis, 2021).

Here, we report, using calculations made with the UK Earth System Model (UKESM1) (Sellar et al., 2020) chemistry–climate model, on the climate and atmospheric composition effects of a range of different atmospheric hydrogen boundary conditions, consistent with a range of different scenarios for a future hydrogen economy. The work presented here is both based on and an extension of the work presented in Warwick et al. (2022). We determine changes in effective radiative forcing (ERF) and present an improved estimate of the hydrogen global warming potential (GWP), accounting for composition changes in both the troposphere and the stratosphere. As part of our GWP calculations, we present a refined methodology for calculating indirect GWPs from steady-state simulations, appropriate for use with both short-lived and longer-lived species.

2 Model and scenarios

2.1 The UK Earth System Model (UKESM1)

The UK Earth System Model (UKESM1) (Sellar et al., 2020), is a state-of-science Earth system model that couples Earth system modules to the HadGEM3-GC3.1 climate model (Kuhlbrodt et al., 2018; Williams et al., 2018). UKESM1 has a horizontal resolution of 1.875 ^∘ in longitude and 1.25 ^∘ in latitude and 85 vertical levels extending from the surface to 85 km. Atmospheric composition changes are calculated using the UK Chemistry and Aerosols (UKCA) atmospheric chemistry module, which includes a coupled stratospheric–tropospheric chemistry scheme (Archibald et al., 2020b) and an interactive two-moment aerosol scheme (Mulcahy et al., 2020). UKCA includes the emissions of nine chemical species: nitric oxide (NO), carbon monoxide (CO), formaldehyde (HCHO), ethane (C₂H₆), propane (C₃H₈), acetaldehyde (MeCHO), acetone (Me₂CO), isoprene (C₅H₈) and methanol (MeOH), while surface mixing ratios of CH₄, N₂O, CFC-11 (CFCl₃), CFC-12 (CF₂Cl₂), CH₃Br, H₂ and COS are prescribed. The photochemical sources and sinks of H₂ and CH₄ are fully interactive, but the use of the lower boundary condition (LBC) fixes the atmospheric burden of H₂ and CH₄. LBCs are widely used in chemistry–climate models as they (1) allow the observed burdens of intermediate lifetime (t > 1 year) species to be imposed, bypassing the need to include emissions, and (2) remove the need for a long spin-up for longer-lived species burdens to reach a steady state. By using an LBC, any chemical feedbacks that would affect the burden of these species are overwritten. However, the effect of feedbacks on the steady-state burden can be computed through quantification of the appropriate chemical feedback factors (see, for example, Heimann et al., 2020, for a discussion on how we have done this for methane).

2.2 Box model

Atmospheric box model simulations were performed using a coupled H₂–CH₄–CO–OH chemical scheme to link lower boundary condition values of H₂ with fugitive H₂ emission rates and obtain CH₄ lower boundary conditions for use in the UKESM1 experiments. In addition, the box model was used for extensive testing of the framework for GWP calculations.

Our box modelling approach is described in more detail in Warwick et al. (2022), is similar to the approach described in Prather (1994) and is used in other studies (Heimann et al., 2020; Nisbet et al., 2020). The model was initialised with realistic values for total methane emissions (585 Tg (CH₄) yr⁻¹), CO (1300 Tg(CO) yr⁻¹) and H₂ (80 Tg(H₂) yr⁻¹) (e.g. Saunois et al., 2020; Pieterse et al., 2013; Zheng et al., 2019) and was found to give methane and H₂ levels in broad agreement with present-day levels – 1865 ppb CH₄, 552 ppb H₂ and 101 ppb CO, as well as a CH₄ lifetime of 9.9 years, which is within the range of observationally derived values.

2.3 Hydrogen economy scenarios

2.3.1 Changes in the abundance of atmospheric H₂ in a hydrogen economy

The increase in abundance of atmospheric H₂ in a future hydrogen economy is not well constrained due to uncertainties in fugitive emissions, which will depend on the ultimate size of the hydrogen economy, H₂ leakage rates and uncertainties in the amount of H₂ undergoing uptake by soils. To guide our model scenarios, we estimate how emissions of H₂ and other species may change in response to a switch from fossil fuel to hydrogen technologies in an illustrative hydrogen economy scenario. In this scenario, approximately 23 % of global energy consumption (about 133 EJ) is supplied by hydrogen (BP, 2020). We, here, assume that in the buildings sector, the transport sector and the power generation sector, 100 %, 50 % and 10 % respectively of the final consumption of fossil fuels switches to hydrogen. The lower percentages for the transport sector and the power generation sector reflect the smaller role hydrogen is assumed to play in these energy sectors due to the existence of low-carbon alternatives such as electric vehicles, wind power and solar power, as well as alternative storage options such as pumped hydro, batteries and compressed air (e.g. Staffell et al., 2019). New hydrogen technologies are assumed to have the same energy efficiency as the fossil fuel technology they are replacing (Staffell et al., 2019), except for in the transport sector, where we assume diesel vehicles and petrol vehicles with an average tank-to-wheel energy efficiency of 30 % will be replaced by vehicles using hydrogen fuel cells with an average efficiency of 50 %. Using a net calorific value of 1 kg H₂ = 33.3 kWh, approximately 860 Tg H₂ yr⁻¹ would be required to provide the energy consumption outlined above.

H₂ leakage rates of the infrastructure associated with a hydrogen economy are uncertain and will depend on many factors (e.g. E4tech, Ltd, 2019). However, all else being equal, H₂ leakage rates of the infrastructure associated with a hydrogen economy are likely to be higher than the leakage rates of natural gas, owing to the small molecule size of hydrogen. A recent study, looking at the US natural gas supply chain, indicated natural gas leaks of around 2.3 % of gross gas production, which is ∼ 60 % higher than the US EPA inventory estimate (Alvarez et al., 2018). H₂ leakage rates of 1 % and 10 % (which represent the range of values used in previous studies e.g. Shultz et al., 2003; Warwick et al., 2004) would lead to additional H₂ emissions of 9 and 96 Tg H₂ yr⁻¹ respectively based on a hydrogen economy supplying 23 % of total present-day energy consumption.

Table 1 List of UKESM1 hydrogen economy scenarios employing different H₂ lower boundary conditions (LBCs). The CH₄ LBC used in the scenarios, not including a CH₄ feedback, is 1835 ppb (2014 year level). In scenarios including a CH₄ feedback (CH4f), the CH₄ LBC is determined using box modelling. Reductions in O₃ precursor emissions (O3pre) are described in Sect. 2.3.2.

We bypass some of the uncertainties in H₂ leakage rates, the ultimate size of the hydrogen economy and soil uptake rates, by performing a range of hydrogen economy scenarios with fixed H₂ lower boundary mixing ratios. We consider scenarios with H₂ lower boundary conditions ranging from 500 to 2000 ppb (i.e. increases relative to present-day mixing ratios of 0 to 1500 ppb), which we believe span many of these uncertainties (see Table 1). For example, based on the size of a hydrogen economy under the conditions outlined above, and assuming the magnitude of the soil sink increases in line with the increase in H₂ mixing ratios (i.e. a constant deposition velocity), box model simulations indicate H₂ lower boundary conditions of 750, 1000 and 1500 ppb represent H₂ leakage rates of about 3 %, 7 % and 13 % respectively (see Sect. 2.2 for further details on the box modelling). However, these H₂ lower boundary conditions would also be consistent with other fractions of total energy supplied by H₂, alternative soil sink responses and H₂ leakage rates. Use of this range of H₂ lower boundary conditions provides a clear signal in the atmospheric response to increased H₂ mixing ratios relative to interannual variability and allows us to explore the linearity of the atmospheric response to increasing H₂. Note that our 2000 ppb H₂ scenario should be considered an extreme end member designed to assess the linearity of the atmospheric response to increasing atmospheric H₂, rather than a projection of potential future atmospheric H₂ levels in a hydrogen economy.

2.3.2 Changes in CH₄, CO, NO_x and VOC emissions in a hydrogen economy

To determine the associated changes in CH₄, CO, NO_x and VOC emissions, further assumptions are required about the new technology employed, in addition to the percentage of different energy sectors switching to H₂. In the building sector, we assume that H₂ will be combusted and NO_x emissions will be limited such that they remain unaltered (despite the potential for higher NO_x emissions from the increased flame temperature), whilst emissions of CH₄, CO and VOCs are eliminated. In the transport sector, we assume a 50 % reduction in emissions of CH₄, CO, NO_x and VOCs based on 50 % of oil-based transport being replaced by hydrogen fuel cells. In the power sector, we assume that stored hydrogen will be used to generate 30 % of the electricity that is currently generated using natural gas via combustion and that NO_x will remain unaltered, whilst emissions of CH₄, CO and VOCs are scaled appropriately. CH₄ emissions associated with energy production are also scaled according to the decreased use of CH₄ where it has been replaced by H₂. The assumptions above lead to reductions in CH₄, CO and NO₂ emissions of 43, 259 and 19 Tg yr⁻¹ (see Warwick et al., 2022, for further details).

Emission reductions for CO, NO_x and VOCs are determined by applying a uniform global-scaling factor to emissions from the specific energy sectors mentioned by Gidden et al. (2019) where hydrogen is assumed to play a role. The emission reductions determined for CH₄ are converted to a change in the CH₄ lower boundary condition for UKESM1 using the box model (see Sect. 2.2). To separate the atmospheric impacts of increasing hydrogen mixing ratios, the feedback on the CH₄ abundance and changes in the emissions of other species emitted by the consumption of fossil fuels, four different sets of hydrogen economy model simulations are performed as described in Table 1.

2.3.3 Model simulations

For the BASE simulation, all boundary conditions are taken from the recently developed Coupled Model Intercomparison Project, Phase 6 (CMIP6), data set and we assume an H₂ LBC of 500 ppb. Climatological boundary conditions for sea-surface temperatures, sea-ice extent and anthropogenic and natural emissions are averaged over the years 2000–2014. In the first set of hydrogen economy scenarios, we consider changes in the atmospheric H₂ abundance only. In these simulations, the methane abundance is fixed at 1835 ppb (the 2014 year level) and does not respond to changes in the methane lifetime. In the second set of hydrogen economy scenarios, we include the methane response to changes in atmospheric H₂ via the impact of H₂ on OH and the methane lifetime. The methane LBCs in these simulations are determined via a series of box model simulations where H₂ is varied (see Sect. 2.2). In the third set of hydrogen economy scenarios, we consider reductions in emissions of the ozone precursors, CO, NO_x and VOCs, from reduced fossil fuel use, whilst the methane abundance is fixed at 1835 ppb. The fourth set of hydrogen economy scenarios includes decreases in the emissions of CH₄, CO, NO_x and VOCs from reduced fossil fuel use, in addition to the feedback of changes in the methane lifetime on the methane abundance. In these scenarios, we assume that all hydrogen will be green hydrogen with no emissions of other species associated with the H₂ production method. In the case that blue or grey hydrogen is used, it is possible methane emissions may increase due to emissions occurring during hydrogen production (e.g. Bertagni et al., 2022). Each simulation was run for 40 years using annually repeating conditions with the final 25 years of each simulation used for analysis and the initial 15 years treated as a spin-up period.

3 Results

3.1 Atmospheric composition impacts following changes in hydrogen and methane

To determine the indirect radiative impact of hydrogen, we need to understand how CH₄, O₃ and H₂O will change, per unit emission of H₂. Those changes will depend not only on how H₂ changes, but also on changes in species emitted from the production and use of fossil fuel. The scenario space is particularly complex and our scenarios should not be regarded as predictions. Instead, here, we try to indicate how the composition of radiatively active species might generally change and which changes are linear in the change in H₂, looking first at the changes driven by H₂ and CH₄.

Figure 1 The response of (a) mass-weighted tropospheric mean OH, (b) CH₄ lifetime with respect to OH, (c) tropospheric O₃ and (d) stratospheric water vapour to increasing H₂ mixing ratios in UKESM1. In panels (a)–(c), the black circle represents the BASE scenario, the blue circles represent scenarios where the CH₄ LBC remains fixed at 2014 levels, the red circles are scenarios where the CH₄ LBC is adjusted to account for the change in CH₄ lifetime, the orange circles are scenarios including changes in emissions of ozone precursors and the green circles include changes in emissions of CH₄ and ozone precursors, as well as adjusted CH₄ LBCs to account for the response of the CH₄ abundance to changes in the CH₄ lifetime. The blue dashed line represents the fit through experiments where the CH₄ LBC remains fixed at 2014 levels and the red dashed line is the fit through experiments which include the response of CH₄ to changing H₂ (and OH). Panel (d) shows the change in the vertical profile of H₂O relative to the BASE scenario where the H₂ LBC is increased to 750, 1000, 1500 and 2000 ppb. Solid lines represent scenarios where CH₄ is fixed at 2014 levels; whereas dashed lines show scenarios where the CH₄ LBC is increased to account for the change in CH₄ lifetime.

Figure 1 summarises the key composition changes seen in a range of different hydrogen lower boundary mixing ratios when emissions of other species from the fossil fuel industry are held constant (blue and red circles) and when they are reduced (orange and green circles). Decreases in OH, the main atmospheric oxidant, are modelled throughout the troposphere with larger decreases in scenarios with higher H₂ abundances (Fig. 1a), consistent with Reaction (R1), being the main atmospheric sink for H₂. The change in OH is linear in H₂ when both the surface CH₄ mixing ratio is held constant (blue circles) and when it is increased to account for the feedback of changes in OH in the methane abundance (red circles). When including this methane feedback, mass-weighted tropospheric mean OH decreases by about 0.90 × 10⁵ molec. cm⁻³ (or about 10 %) for every 1000 ppb (1 ppm) increase in H₂ (Fig. 1a, red circles). When tropospheric mean OH is weighted by reaction with CH₄, OH decreases by about 1.8 × 10⁵ molec. cm⁻³ for every 1 ppm increase in H₂. These modelled changes in OH go on to cause a cascade of further composition changes by altering the methane lifetime, production and destruction rates of tropospheric ozone, as well as aerosol nucleation rates and cloud condensation nuclei (CCN) concentrations.

Figure 1b shows the CH₄ lifetime with respect to loss by reaction with tropospheric OH as a function of atmospheric H₂ and tropospheric OH. In the BASE simulation, the modelled CH₄ lifetime with respect to loss by reaction with tropospheric OH is 8.5 years. Adjusting to include losses in the stratosphere (120-year lifetime), but excluding soil loss (160-year lifetime), gives a chemical CH₄ lifetime of 7.9 years. This value is lower than the range given for the total chemical methane lifetime in IPCC AR6 (9.7 ± 1.1 years, Szopa et al., 2021) and a present-day total chemical multimodel estimate of 8.4 ± 0.3 years in Stevenson et al. (2020). As atmospheric H₂ increases but all other factors remain constant, there is a linear increase in the CH₄ lifetime. The methane lifetime increases by 0.64 years for every 1000 ppb (1 ppm) increase in H₂ when the methane abundance is held constant (blue circles). However, the methane lifetime increases by 0.96 years for every 1000 ppb increase in H₂ when the methane abundance is increased in response to the modelled changes in its lifetime (red circles).

The oxidation of H₂ in the troposphere can also affect the levels of ozone (O₃) as the HO₂ produced through Reaction (R1) can catalyse the interconversion of NO_x which drives tropospheric ozone production (e.g. as reviewed in Archibald et al., 2020a),

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

This is followed by photolysis of NO₂ to form ozone,

$$\begin{array}{}\text{(R4)}& {\mathrm{NO}}_{\mathrm{2}}+h\mathit{\nu }\to \mathrm{NO}+\mathrm{O}(+{\mathrm{O}}_{\mathrm{2}}={\mathrm{O}}_{\mathrm{3}}).\end{array}$$

In addition, changes in OH and HO₂ produced via Reaction (R1) influence the destruction of tropospheric ozone by altering the flux through the reaction of O₃ with both of these species. Overall, the tropospheric ozone burden increases approximately linearly with atmospheric H₂ increase, driven primarily through increases in the reaction HO₂ + NO (Reaction R3), which increases by ∼ 7 % when atmospheric H₂ increases to 2000 ppb (Fig. 1c). Changes in the ozone budget are documented in Table S2 in the Supplement.

Increased H₂ mixing ratios at the surface can also influence stratospheric composition. Increases in H₂ and CH₄ abundance in the troposphere result in an increased flux of these species to the stratosphere, where their oxidation leads to the production of water vapour. Figure 1d shows the increase in the stratospheric water vapour profile relative to the BASE scenario when tropospheric H₂ increases from 500 to 750, 1000, 1500 and 2000 ppb, both where the tropospheric methane abundance remains unchanged and where the methane LBC is adjusted in response to modelled changes in the methane lifetime. When H₂ increases to 2000 ppb but CH₄ is held constant, stratospheric H₂O increases by up to 25 % (> 1 ppm). Note that in our simulations the maximum increase in stratospheric water vapour (occurring at ∼ 55 km) is slightly less than 1 ppm per 1000 ppb increase in surface H₂. If the tropospheric methane abundance is increased in line with modelled changes in the CH₄ lifetime, larger increases in stratospheric H₂O are simulated.

3.2 Atmospheric composition impacts when hydrogen, methane and ozone precursor emissions all change

The impacts discussed above arise from changes in atmospheric H₂ abundance and subsequent changes in methane abundance, arising from the feedback of OH on the methane lifetime. Expected reductions in ozone precursor emissions of CO, NO_x (which may also increase locally), VOCs and CH₄, following a shift to a global hydrogen economy, also have the potential to further impact atmospheric composition. The green circles and the orange circles in Fig. 1a–c show the impact of including these other emission changes with and without the feedback on the methane abundance included respectively. The relationship remains linear for the air-mass-weighted tropospheric mean OH, methane lifetime and tropospheric ozone burden plotted as a function of H₂ LBC. In all cases, the gradients remain similar to the corresponding scenarios not including the ozone precursor emission reductions (note the ozone precursor emission changes are the same across all scenarios).

The linear relationships in Fig. 1 show allow the global methane and ozone changes in a range of different hydrogen economy scenarios to be estimated, giving responses to different realisations of a future hydrogen economy. They can also be employed to obtain a GWP for hydrogen (see Sect. 3.4).

Figure 2 Shows the effect of increases in atmospheric H₂ on (a) tropospheric ozone column for the BASE scenario (DU); (b) the change in tropospheric ozone column relative to BASE when H₂ is increased to 2000 ppb, including the response of the methane abundance to changes in the methane lifetime (DU); (c–f) changes in tropospheric ozone column relative to BASE when H₂ mixing ratios are increased, but all other factors including the methane abundance remain fixed; and (g–j) changes in tropospheric ozone column relative to BASE for scenarios including emission reductions in ozone precursors and methane. The numbers above each figure give the change in global mean tropospheric column in DU relative to BASE. Only differences that are statistically significant at the 95 % confidence level are presented.

Some detail on geographical changes in tropospheric ozone in a hydrogen economy is provided in Fig. 2. As the photochemical production of ozone is driven by oxidation of CO, CH₄ and VOCs in the presence of NO_x, the net regional change in ozone depends not only on how H₂ changes, but also on changes in emissions of these species. This sensitivity is demonstrated in Fig. 2.

When H₂ increases, so does tropospheric ozone (Fig. 2c to f; see also Fig. 1c). Similarly, when atmospheric methane increases, the effect is to enhance the tropospheric ozone column increase (compare Fig. 2b with f). In contrast, reductions in ozone precursor and CH₄ emissions would avoid the significant O₃ increases seen in Fig. 2b (see Fig. 2h and j). For example, if there is no H₂ leakage, reductions in the other emissions (Fig. 2g and h) lead, as expected, to reductions in tropospheric ozone. When surface H₂ reaches 2000 ppb, reductions in the other emissions still lead to modest decreases in ozone (Fig. 2i, j). For example, in the UKESM1 simulation, which assumes a large H₂ leakage, with surface mixing ratios of H₂ reaching 2000 ppb and emission reductions in CH₄, CO, NO_x and VOCs, the global mean tropospheric column ozone response is found to be small (−0.1 DU) due to the competing effects outlined above (Fig. 2j). The tropospheric ozone response to a global shift towards a hydrogen economy is therefore strongly influenced by the amount of hydrogen added to the atmosphere through leaks and the co-benefit reductions achieved in CO, NO_x, VOC and CH₄ emissions.

As shown above, the changes in H₂ and methane affect OH which can in turn also influence aerosols in our model simulations. We simulate a linear decrease in both tropospheric air-mass-weighted OH (Fig. 1a) and the OH+SO₂ reaction flux, for increasing H₂ in the atmosphere. The reduction in the flux through the OH + SO₂ reaction results in a shift from the oxidation of SO₂ in the gas phase (which contributes towards new particle formation) towards the oxidation of SO₂ in the condensed phase (which grows existing particles). Such a shift leads to a reduction in the aerosol number and an increase in aerosol size with impacts on clouds and radiative forcing (see O'Connor et al., 2022, for further details). The modelled decrease in OH following the increase in H₂ is coincident with an increase in HO₂ (R1 leads to the production of H atoms, which in the troposphere near instantaneously form HO₂); thus, H₂ acts to increase the HO₂ : OH ratio, similar to CO. We find that the HO₂ : OH ratio linearly increases with increasing H₂ in the atmosphere (from 145 at 500 ppb H₂ to 175 at 2000 ppb H₂), and, as a result, the H₂O₂ concentration increases, further enhancing in-cloud oxidation of SO₂. Whilst our simulations demonstrate the potential importance of aerosol feedbacks relating to changes in OH, the inclusion of other aerosols, e.g. nitrates, may influence our results and more studies, involving multiple models, are required to constrain the uncertainties involved.

Figure 3 The effect of the simulated changes in atmospheric composition on effective radiative forcing (ERF). Panel (a) shows the ERF as a function of H₂ and other experimental conditions and panel (b) shows the ERF as a function of CH₄ and other experimental conditions. Just as for Fig. 1a–c, the black circle represents the BASE scenario, the blue circles represent scenarios where the CH₄ LBC remains fixed at 2014 levels, the red circles are scenarios where the CH₄ LBC is adjusted to account for the change in CH₄ lifetime, the orange circles are scenarios including changes in emissions of ozone precursors, and the green circles include changes in emissions of CH₄ and ozone precursors, as well as adjusted CH₄ LBCs to account for the response of the CH₄ abundance to changes in the CH₄ lifetime. Note that none of the calculated ERFs include anticipated reductions in CO₂.

Table 2 ERFs and their clear-sky and cloud radiative effect components in the shortwave (SW) and longwave (LW) calculated using the method outlined in Ghan (2013, as used in O'Connor et al., 2021, 2022), based on the last 30 years of hydrogen simulations. All units are in watts per square metre (W m⁻²).

3.3 Effective radiative forcing

Figure 3 shows the effective radiative forcing (ERF) for the scenarios in Table 1. Our ERF values are based on the chemical changes in composition resulting from increases in atmospheric H₂, in addition to changes in emissions of CH₄, CO, NO_x and VOCs in the scenarios where these are included. The calculated ERFs do not include the large climate benefit expected from anticipated reductions in CO₂ emissions. We determined the ERF following Eq. (8) from O'Connor et al. (2021). The ERF can also be decomposed into a clear-sky and a cloud radiative forcing, with the latter calculated from a change in the cloud radiative effect (CRE) from “clean” radiation calls that exclude aerosol–radiation interactions (Ghan, 2013). The clear-sky forcing, in principle, may include a contribution from aerosol–radiation interactions, but the magnitude of this contribution was found to be less than 0.02 W m⁻². The time-slice experiments using an atmosphere-only model, i.e. using decoupled sea surface temperature (SST) and sea-ice coverage, permit the determination of the top-of-atmosphere radiative forcing including rapid adjustments to cloud and water vapour. The ERF was calculated relative to the 2014 base case. The ERF varies from a cooling (negative ERF) to a warming (positive ERF) tendency, depending on the scenario. The ERF increases with H₂ over the range 500–2000 ppb. Table 2 shows that these increases are due to both clear-sky and cloud radiative forcing. The clear-sky forcing increases with increasing H₂, presumably due to the increased O₃, while the CRE can be ascribed to small decreases in cloud albedo. Figure 4 shows that increasing H₂ from 500 to 2000 ppb leads to changes in cloud droplet number concentration (CDNC) across the globe. The decreased CDNC leads, for the same amount of water vapour, to larger cloud droplets and lower cloud albedo. This leads to a CRE of 0.05 W m⁻². In scenarios where CH₄ is increased, the further suppression of OH leads to a stronger CRE whilst the radiative forcing from the CRE in our hydrogen economy scenarios increases the total ERFs by approximately 50 %, thus forming a significant component of indirect forcing from H₂. We do not include forcing from clouds in our GWP calculations. Such couplings between oxidants, aerosols and clouds are relatively unique in the UKESM1 compared with other CMIP6 models, and more studies are required to constrain the uncertainties involved.

Figure 4 Cloud droplet number concentration (CDNC) in the UKESM1 model at 1 km altitude. Panel (a) shows CDNC in the BASE (500 ppb H₂) scenario, panel (b) the change in CDNC when H₂ is increased to 2000 ppb (2000H2 scenario).

3.4 Global warming potential

A GWP quantifying the indirect radiative forcings associated with H₂ has, so far, only been determined in a limited number of studies previously, including studies by Derwent et al. (2001, 2020, an update on Derwent et al., 2006), Field and Derwent (2021), Hauglustaine et al. (2022), and Sand et al. (2023). The first two studies considered only the influence of H₂ on tropospheric composition, whereas Hauglustaine et al. (2022) and Sand et al. (2023) included both the tropospheric and stratospheric response. In addition to these studies, Ocko and Hamburg (2022) extended the work presented in Warwick et al. (2022), to calculate hydrogen's GWP for a range of time horizons.

We determine a GWP for H₂ based on composition changes in both the troposphere and stratosphere in our UKESM1 simulations, using radiative forcing scaling factors from Myhre et al. (2013), Forster et al. (2021) and Eqs. (1) to (4) presented in Sect. 3.4.1 (see also Warwick et al., 2022). Our new estimate of the H₂ GWP considers changes in methane and ozone in the troposphere and changes in stratospheric water vapour and stratospheric ozone. Note that we have not used the model-derived ERFs from Sect. 3.3 in the GWP calculations below. The model-derived ERFs (for H₂ leakage only, including the methane response) are larger than the radiative forcing due to hydrogen determined from modelled changes in chemical composition using radiative forcing scaling factors from Forster et al. (2021) (see Sect. 3.4.2). Much of this difference arises due to the large contribution from clouds to the model-based ERFs (see Table 2) that is not accounted for in the radiative forcing scaling factors. For our GWP calculations, we determine the radiative forcing using the modelled changes in composition along with the radiative forcing scaling factors rather than the ERFs. Our reason for this is two-fold. Firstly, it ensures the large uncertainty associated with the aerosol–cloud component of the ERFs (Forster et al., 2021) is not propagated into the GWP calculations. Secondly, significant further work would be required to attribute the calculated total ERFs to individual species so that the appropriate lifetimes could be applied to the different components of the forcing that contribute to the GWP. However, we note that if the model ERFs were to be used in the GWP calculations, the GWPs presented in the next section would be larger.

3.4.1 Methodology

For our GWP calculation, we derive a more universal version of the approach presented by Fuglestvedt et al. (2010) for calculating GWPs for species whose emissions result in indirect radiative forcings, which was also implicitly adopted for global temperature change potentials in the IPCC Sixth Assessment Report (Forster et al., 2021). The method of Fuglestvedt et al. (2010) assumes that the time evolution of the radiative forcing depends only on the lifetime of the species causing the forcing and not the lifetime of the species being emitted. In their method, the time evolution of the radiative forcing during a 1-year constant emission is described by an approach to a steady state with a time constant according to the lifetime of the species causing the forcing. The peak forcing is assumed to occur directly at the end of the 1-year emission and to subsequently decay with the same lifetime as in the growth phase. Although this assumption may be a good approximation for short-lived emission species such as NO_x, it may not be for species with longer lifetimes such as H₂. Our extension of their method (see Warwick et al., 2022) accounts for both the lifetime of the emitted species and the lifetime of the species causing the forcing, as well as allowing the time period of the constant emission to be varied. The use of a 1-year constant emission rather than an emission pulse has the benefit of removing the time-of-year dependence for when the emission pulse occurs.

In our method for deriving indirect GWPs from steady-state simulations, we include the indirect forcing due to perturbations in radiatively active species arising (a) during the constant emission (AGWP1, Eq. 1), (b) from the decay of the perturbation subsequent to the end of the constant emission (AGWP2, Eq. 2) and (c) from perturbations arising from the emission species remaining in the atmosphere after the end of the constant emission (AGWP, Eq. 3). Derivations for Eqs. (1) to (3) are presented in the Supplement.

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}\mathrm{AGWP}\mathrm{1}& ={\displaystyle \frac{R{a}_M{\mathit{\alpha }}_M{\mathit{\alpha }}_{\mathrm{H}}C}{t_p}}\left(t_p-{\mathit{\alpha }}_M\left(\mathrm{1}-\exp \left({\displaystyle \frac{-t_p}{{\mathit{\alpha }}_M}}\right)\right)\right.\\ & -\left({\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{H}}}{{\mathit{\alpha }}_{\mathrm{H}}-{\mathit{\alpha }}_M}}\right)\left({\mathit{\alpha }}_{\mathrm{H}}\left(\mathrm{1}-\exp \left({\displaystyle \frac{-t_p}{{\mathit{\alpha }}_{\mathrm{H}}}}\right)\right)\right.\\ & \left.\left.-{\mathit{\alpha }}_M\left(\mathrm{1}-\exp \left({\displaystyle \frac{-t_p}{{\mathit{\alpha }}_M}}\right)\right)\right)\right),\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}\mathrm{AGWP}\mathrm{2}& ={\displaystyle \frac{\left(R{a}_M{\mathit{\alpha }}_{\mathrm{H}}^{\mathrm{2}}{\mathit{\alpha }}_{M}C\left(\mathrm{1}-\exp \left(\frac{-t_p}{{\mathit{\alpha }}_{\mathrm{H}}}\right)\right)\right)}{t_p\left({\mathit{\alpha }}_{\mathrm{H}}-{\mathit{\alpha }}_M\right)}}\\ & \left({\mathit{\alpha }}_{\mathrm{H}}\left(\mathrm{1}-\exp \left(-{\displaystyle \frac{\left(\mathrm{Hz}-t_p\right)}{{\mathit{\alpha }}_{\mathrm{H}}}}\right)\right)\right.\\ & \left.-{\mathit{\alpha }}_M\left(\mathrm{1}-\exp \left(-{\displaystyle \frac{\left(\mathrm{Hz}-t_p\right)}{{\mathit{\alpha }}_M}}\right)\right)\right),\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}\mathrm{AGWP}\mathrm{3}& ={\displaystyle \frac{R{a}_M{\mathit{\alpha }}_M^{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{H}}C}{t_p}}\left(\left(\mathrm{1}-\exp \left(-{\displaystyle \frac{t_p}{{\mathit{\alpha }}_M}}\right)\right)\right.\\ & \left.-\left({\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{H}}}{{\mathit{\alpha }}_{\mathrm{H}}-{\mathit{\alpha }}_M}}\right)\left(\exp \left(-{\displaystyle \frac{t_p}{{\mathit{\alpha }}_{\mathrm{H}}}}\right)-\exp \left(-{\displaystyle \frac{t_p}{{\mathit{\alpha }}_M}}\right)\right)\right)\\ & \left(\mathrm{1}-\exp \left(-{\displaystyle \frac{\left(\mathrm{Hz}-t_p\right)}{{\mathit{\alpha }}_M}}\right)\right),\end{array}\end{array}$$

where AGWP is the absolute global warming potential (W m⁻² kg⁻¹ yr); M is the species resulting in the indirect radiative forcing, CH₄, O₃ and H₂O (ppb or DU for tropospheric O₃); R_M is the radiative forcing efficiency for M (W m⁻² ppb⁻¹ or W m⁻² DU⁻¹ for tropospheric O₃); a_M is the production rate of M (ppb yr⁻¹) per ppb H₂ change at the steady state (yr⁻¹); α_M is the atmospheric lifetime of M (years); α_H is the atmospheric H₂ lifetime (combined chemical and deposition lifetime) (years); Hz is the time horizon considered (years); C is the conversion factor for converting H₂ mixing ratio (ppb) into H₂ mass (kg) and t_p is the length of step emission (years).

Equations (1) to (3) are applied separately to methane, non-methane-induced ozone and stratospheric water vapour to calculate a net GWP for H₂. The associated indirect GWPs for each radiative species perturbed by hydrogen are then given by

$$\begin{array}{}\text{(4)}& {\mathrm{GWP}}_{{\mathrm{H}}_{\mathrm{2}}}={\displaystyle \frac{{\mathrm{AGWP}}_{{\mathrm{H}}_{\mathrm{2}}}}{{\mathrm{AGWP}}_{{\mathrm{CO}}_{\mathrm{2}}}}}.\end{array}$$

A comparison of the time evolution of the methane perturbation resulting from a change in H₂ mixing ratio based on the standard (Fuglestvedt et al. 2010) equations and Fig. 5 shows our updated equations. Accounting for the atmospheric lifetime of hydrogen in the updated equations results in a slower predicted rate of increase in methane during the 1-year constant emission relative to the standard (Fuglestvedt et al. 2010) equations, a delayed methane peak (peaking at 5 years rather than 1 year, using the parameters specified in the caption) and a slower rate of methane decrease subsequent to that peak, as methane is still being impacted by the remaining H₂. This difference occurs because the equations of Fuglestvedt et al. (2010) assume that the emitted species instantaneously reaches a new steady-state atmospheric concentration at the start of the 1-year constant emission. In addition, their assumption that the decay following the peak forcing is controlled only by the lifetime of the radiatively active species neglects any subsequent perturbations to atmospheric chemistry as a result of the emitted species still present in the atmosphere after the end of the 1-year constant emission. A comparison of the transient methane growth due to a H₂ perturbation, as determined by Eqs. (1) to (3) in Fig. 5 (grey line), with model studies where methane has been free to respond to a H₂ perturbation (including the box model described in Sect. 2.2 and the model studies of Derwent et al., 2020, and Bertagni et al., 2022) shows that the response represented by Eqs. (1) to (3) is similar to that found in transient model simulations following a H₂ pulse emission.

Figure 5 The time evolution of excess methane during and after a 1-year step hydrogen emission (based on a steady-state H₂ excess of +1500 ppb, an H₂ lifetime of 1.96 years and a methane perturbation lifetime of 11.8 years) as predicted by the equations of Fuglestvedt et al. (2010) (yellow) and our updated equations. For our updated equations, the excess methane is split into contributions from methane generated during the step emission (blue), methane generated from excess hydrogen remaining in the atmosphere after the 1-year step emission (orange). The total methane excess in shown in grey.

Equations (1) to (3) are obtained using a framework assuming a 1-year constant emission rather than an emission pulse. This has the benefit of removing the time-of-year dependence of the atmospheric perturbations arising from the emission pulse. However, once obtained, the chemical parameters from our UKESM1 simulations can also be applied to the pure pulse approach with little error. Equations (1) to (3) can be simplified by using the t_p→0 limit of Eq. (2), which is given in Eq. (5). Based on our obtained parameters, using Eq. (5) rather than Eqs. (1) to (3) leads to a difference of 0.1 % for our hydrogen GWP(100) and 1.5 % for our hydrogen GWP(20). In contrast, using the original method of Fuglestvedt et al. (2010) leads to a 0.1 % difference in our hydrogen GWP(100) but a 5 % difference in our hydrogen GWP(20) relative to Eqs. (1) to (3).

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\mathrm{AGWP}& ={\displaystyle \frac{R{a}_M{\mathit{\alpha }}_{\mathrm{H}}{\mathit{\alpha }}_{M}C}{\left({\mathit{\alpha }}_{\mathrm{H}}-{\mathit{\alpha }}_M\right)}}\left({\mathit{\alpha }}_{\mathrm{H}}\left(\mathrm{1}-\exp \left(-{\displaystyle \frac{\mathrm{Hz}}{{\mathit{\alpha }}_{\mathrm{H}}}}\right)\right)\right.\\ & \left.-{\mathit{\alpha }}_M\left(\mathrm{1}-\exp \left(-{\displaystyle \frac{\mathrm{Hz}}{{\mathit{\alpha }}_M}}\right)\right)\right)\end{array}\end{array}$$

3.4.2 Calculation of a hydrogen GWP

In our calculations, AGWP${}_{{\mathrm{CO}}_{\mathrm{2}}}$(100) and AGWP${}_{{\mathrm{CO}}_{\mathrm{2}}}$(20) are taken to be 9.17 × 10⁻¹⁴ (± 26 %) and 2.49 × 10⁻¹⁴ (± 18 %) W m⁻² kg⁻¹ yr (Myhre et al., 2013) respectively. AGWP${}_{{\mathrm{H}}_{\mathrm{2}}}$ is taken as the sum of AGWP1–3. The values for AGWP${}_{{\mathrm{CO}}_{\mathrm{2}}}$(100) and AGWP${}_{{\mathrm{CO}}_{\mathrm{2}}}$(20) change slightly to 8.95 × 10⁻¹⁴ (± 26 %) and 2.43 × 10⁻¹⁴ (± 18 %) W m⁻² kg⁻¹ yr respectively in Forster et al. (2021). Table 3 shows the impact of this on our H₂ GWP calculation, along with other IPCC AR6 (Forster et al., 2021) updates. The length of the step emission, t_p, is 1 year, Hz is taken to be 20 or 100 years, and a_H is 1.96 years (with an uncertainty range from 1.4 to 2.2 years). The conversion factor, C, for converting H₂ mixing ratio into H₂ mass based on the UKESM1 data is 3.52 × 10⁻⁹ ppb kg⁻¹. Values for R_M, a_M and α_M are dependent on the species causing the indirect forcing (methane, ozone and stratospheric water vapour).

Table 3 A comparison of the 100- and 20-year time horizon H₂ GWPs (GWP(100) and GWP(20) respectively) determined using radiative efficiencies and the AGWP for CO₂ from IPCC AR5 (Myhre et al., 2013) and from IPCC AR6 (Forster et al., 2021). The uncertainties quoted account for uncertainties in the H₂ lifetime with respect to soil uptake, radiative forcing scaling factors and the AGWP${}_{{\mathrm{CO}}_{\mathrm{2}}}$ only.

For methane, $R_{{\mathrm{CH}}_{\mathrm{4}}}$ is taken to be 0.000363 (± 10 %) W m⁻² (ppb CH₄)⁻¹ (Myhre et al., 2013) or 0.000389 W m⁻² (ppb CH₄)⁻¹ (Forster et al., 2021). The methane perturbation lifetime, ${\mathit{\alpha }}_{{\mathrm{CH}}_{\mathrm{4}}}$, in UKESM1 is determined to be 11.8 years in this study, which falls within the range given by IPCC AR5 (12.4 ± 1.4 years, Myhre et al., 2013). The value for $a_{{\mathrm{CH}}_{\mathrm{4}}}$, an additional methane production rate term per unit increase in H₂, represents the positive CH₄ tendency resulting from an increase in H₂ and corresponding decrease in OH and is defined as

$$\begin{array}{}\text{(6)}& a_{{\mathrm{CH}}_{\mathrm{4}}}={\mathrm{CH}}_{\mathrm{4}}{k}_{{\mathrm{CH}}_{\mathrm{4}}+\mathrm{OH}}{\displaystyle \frac{\mathrm{d}\mathrm{OH}}{\mathrm{d}{\mathrm{H}}_{\mathrm{2}}}},\end{array}$$

where CH₄ is the methane lower boundary condition (1835 ppb in UKCA), $k_{{\mathrm{CH}}_{\mathrm{4}}+\mathrm{OH}}$ is the global mean air-mass-weighted rate constant for the reaction of CH₄ with OH and $\mathrm{d}\mathrm{OH}/\mathrm{d}{\mathrm{H}}_{\mathrm{2}}$ is the rate of change in air-mass-weighted tropospheric mean OH with respect to H₂. Values for these parameters derived from UKESM1 simulations are as follows: k(CH₄+OH) = 2.8 × 10⁻¹⁵ cm³ molec.⁻¹ s⁻¹ and $\mathrm{d}\mathrm{OH}/\mathrm{d}{\mathrm{H}}_{\mathrm{2}}$ = 3.68 × 10⁻⁹ (dimensionless). k(CH₄+OH) is determined from the gradient of methane lifetime against 1 $/$ (global air-mass-weighted mean OH) and $\mathrm{d}\mathrm{OH}/\mathrm{d}{\mathrm{H}}_{\mathrm{2}}$ from the gradient of the mass-weighted tropospheric mean OH concentration against H₂ surface mixing ratio for the UKESM1 simulations including the feedback on the methane lifetime via an adjusted CH₄ lower boundary condition determined by the box model (see Fig. 1a). Note that if we weight OH by reaction with CH₄ instead of air mass, the value determined for k(CH₄+OH) halves (1.4 × 10⁻¹⁵ cm³ molec.⁻¹ s⁻¹), the value for $\mathrm{d}\mathrm{OH}/\mathrm{d}{\mathrm{H}}_{\mathrm{2}}$ doubles (7.36 × 10⁻⁹ (dimensionless)) and our value for $a_{{\mathrm{CH}}_{\mathrm{4}}}$ remains unaltered. It is also possible to derive $\mathrm{d}\mathrm{OH}/\mathrm{d}{\mathrm{H}}_{\mathrm{2}}$ using the gradient of the mass-weighted tropospheric mean OH concentration against H₂ surface mixing ratio from the UKESM1 simulations where only the H₂ LBC is varied (whilst the CH₄ LBC is held constant), in combination with the equation below from Stevenson et al. (2013), and the methane feedback factor of 1.39 derived for UKESM1,

$$\begin{array}{}\text{(7)}& {\mathrm{CH}}_{\mathrm{4}\mathrm{new}}={\mathrm{CH}}_{\mathrm{4}\mathrm{base}}{\left({\displaystyle \frac{{\mathit{\tau }}_{\mathrm{new}}}{{\mathit{\tau }}_{\mathrm{base}}}}\right)}^f,\end{array}$$

where f is the methane feedback factor, CH_4 base and τ_base are the lower boundary condition and methane lifetime in the base scenario, and CH_4 new and τ_new are the lower boundary condition and methane lifetime in the perturbed H₂ scenario. The difference between the $\mathrm{d}\mathrm{OH}/\mathrm{d}{\mathrm{H}}_{\mathrm{2}}$ values obtained using the two methods above is approximately 3 %. This, in turn, leads to a 3 % difference in the methane-induced hydrogen GWP terms and a difference of less than 1 % in the total hydrogen GWP when accounting for all forcings. This difference is negligible in comparison to other uncertainties (see Table 3); therefore our GWP calculation remains independent of the box model results.

The indirect forcing from methane-induced changes in tropospheric ozone and stratospheric water vapour are determined by multiplying the calculated direct methane forcing ($R_{{\mathrm{CH}}_{\mathrm{4}}}$) by scaling factors of 0.5 ± 55 % and 0.15 ± 70 % respectively (Myhre et al., 2013). Equivalent scaling factors, from Forster et al. (2021), are 0.37 ± 50 % and 0.106 ± 100 %.

In addition to the methane-induced indirect radiative forcing from ozone and water vapour outlined above, we also consider non-methane-induced radiative forcing from H₂ as a source of tropospheric ozone and stratospheric water vapour. For stratospheric water vapour, the e-folding lifetime of the water vapour perturbation following a change in H₂, ${\mathit{\alpha }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$, is derived from model spin-up data following a change in the H₂ lower boundary condition and is determined to be 8 years in the lower stratosphere. Values for both $R_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ and $a_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ are determined by the steady-state change in the stratospheric water profile between our UKESM1 BASE simulation and the simulation where the H₂ lower boundary condition is increased from 500 to 2000 ppb. A comparison of the vertical profile of stratospheric water vapour changes resulting from H₂ changes in our study (Fig. 1d), with that resulting from methane changes between 1950 and 2000 in Fig. 2b of Myhre et al. (2007), shows very similar H₂O changes at all altitudes. We therefore adopt the radiative forcing determined by Myhre et al. (2007) of 0.05 W m⁻² with an estimated uncertainty of ±20 % for stratospheric water vapour changes relative to base in our 2000H₂ UKESM1 scenario. This water vapour forcing is associated with a change of ∼ 500 ppb in stratospheric water vapour at 30 km altitude, giving a value for $R_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ of 1 × 10⁻⁴ W m⁻² (ppb H₂O)⁻¹ when considering changes at 30 km. We assume that the entire stratospheric H₂O profile will scale proportionally with the H₂ lower boundary condition (following results from our UKESM1 simulations using H₂ LBCs of 750, 1000, 1500 and 2000 ppb H₂), so the choice of 30 km for determining $R_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ is arbitrary. The production rate of water vapour per ppb change in H₂, $a_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$, can then be defined as the change in steady-state H₂O at 30 km (500 ppb) divided by the lifetime of 8 years and the change in H₂ mixing ratio of 1500 ppb.

For changes in the global mean tropospheric ozone column, the value for $R_{{\mathrm{O}}_{\mathrm{3}}}$ is 0.042 (0.037 to 0.047) W m⁻² DU⁻¹ and is taken directly from Myhre et al. (2013). No updated value for $R_{{\mathrm{O}}_{\mathrm{3}}}$ was presented by Forster et al. (2021), although a similar value of 0.039 W m⁻² DU⁻¹ can be obtained based on the estimated change in tropospheric ozone of 109 ± 25 Tg between 1850 and 2005 and a stratospheric-temperature-adjusted radiative forcing due to tropospheric ozone of 0.39 ± 0.07 W m⁻² between 1850 and 2010. The e-folding lifetime associated with the ozone perturbation, ${\mathit{\alpha }}_{{\mathrm{O}}_{\mathrm{3}}}$, is taken to be 0.07 years, calculated as the tropospheric O₃ burden divided by the loss in our UKESM1 base simulation. This is shorter than the lifetime used for the short-lived ozone perturbation by Fuglestvedt et al. (2010) of 0.267 years. However, our results are insensitive to the difference between these values and small variations around them. The value for $a_{{\mathrm{O}}_{\mathrm{3}}}$ is calculated as the change in the steady-state global mean tropospheric ozone column in our 2000H₂ simulation relative to base (0.59 DU) divided by the lifetime of 0.07 years and the change in H₂ mixing ratio of 1500 ppb. A list of all the input parameters, with the values used in our AGWP calculations, is presented in Table S1. Changes in the H₂ mixing ratio can also influence stratospheric ozone. In our UKESM1 simulations, the change in global mean stratospheric ozone column between our base and 2000H₂_CH4f (increased H₂, including the chemical feedback on the methane lifetime) scenarios was negligible at +0.10 DU, relative to a base value of 283.78 DU and interannual variations of up to several DU. Based on a value of R for stratospheric ozone of 0.0054 W m⁻² DU⁻¹ (Schwarzkopf and Ramaswamy, 1993) and a range of e-folding lifetimes for the stratospheric ozone perturbation of between 1 and 10 years, the results obtained are a GWP(100) of less than 0.03 using Eqs. (1) to (4). We consider this to be below the uncertainty involved in the calculation and conclude that changes in stratospheric ozone for our range of scenarios do not significantly contribute to the hydrogen GWP.

For a 100-year time horizon we determine a H₂ GWP of 12 ± 6 (see Table 3). This value is larger than previously published studies, which do not account for changes in stratospheric composition (Derwent et al., 2020, 5 ± 1, and Field and Derwent, 2021, 3.3 ± 1.4), but similar to the values of Hauglustaine et al. (2022) (12.8 ± 5.2) and the multimodel study of Sand et al. (2023) (11.6 ± 2.8, 1 standard deviation), which include a tropospheric and stratospheric response. Approximately two-thirds of our GWP arise from the influence of H₂ on methane and ozone distributions in the troposphere, and one-third arises from the influence of H₂ on stratospheric water vapour. Results from our study can also be compared with that of Paulot et al. (2021), who estimated the indirect radiative forcing at a steady state due to an increase in H₂. Their study accounted for changes in methane, tropospheric ozone and stratospheric water vapour and found an indirect radiative forcing arising from a H₂ of 0.13 mW m⁻² (ppb H₂)⁻¹. Using the radiative forcing scaling factors, we obtain a larger value of 0.18 ± 0.03 mW m⁻² (ppb H₂)⁻¹ (the same value is obtained using both the Myhre et al., 2013, and Forster et al., 2021, radiative efficiency scaling factors), which may partly be explained by the larger methane feedback factor in UKESM1 (1.39) relative to the chemistry–climate model used in (Paulot et al., 2021).

Uncertainties in our calculation are based on uncertainties in the hydrogen lifetime with respect to soil uptake, uncertainties in radiative forcing scaling factors and the AGWP for CO₂. We do not account for uncertainties in chemical lifetimes and other UKESM1-derived quantities which are likely to vary when using different atmospheric models with different chemistry schemes. In addition, we also note that the sensitivity of ERF in UKESM1 to changes in methane is approximately 65 % larger than that indicated by the radiative scaling factors used in the empirical relationships applied in the GWP calculation (O'Connor et al., 2022). This suggests a further source of uncertainty than accounted for by the radiative forcing scaling factor uncertainties outlined by Myhre et al. (2013) and Forster et al. (2021), associated with the response of aerosols, clouds, ozone and water vapour to changing methane levels.

4 Discussion and conclusions

The leakage of hydrogen associated with a hydrogen economy will result in indirect global warming, offsetting greenhouse gas emission reductions made as a result of a switch from fossil fuel to hydrogen. We have presented a methodology for calculating the indirect GWP of H₂, with principal contributions coming from changes in methane, tropospheric ozone and stratospheric water vapour. The GWP depends on the lifetime of H₂ in the atmosphere, the perturbation lifetimes of the radiatively active species (e.g. CH₄ and O₃) and their effective production rates per unit change in H₂. We present values for these based on calculations with our chemistry–climate model, UKESM1 and a box model. The 100-year GWP for H₂, based on the above, is 12 ± 6, a value somewhat higher than some previous estimates but similar to recent values from Hauglustaine et al. (2022) and Sand et al. (2023). Whilst our GWP uncertainty accounts for uncertainties associated with the size of the soil sink and radiative forcing scaling factors, it does not include uncertainties in model-derived parameters (e.g. chemical and perturbation lifetimes or the chemical production rates per unit change in H₂). In particular, uncertainties associated with the stratospheric response may be significant, as well as the response of aerosols and clouds to changing methane levels. Further simulations using different Earth system models and different chemistry schemes would be required to fully assess the impact of uncertainties in these parameters on the hydrogen GWP.

Determining a GWP for hydrogen allows a change in hydrogen emissions to be compared to an equivalent change in carbon dioxide emissions in terms of time-integrated radiative forcing. This increase in equivalent carbon dioxide emissions can be compared with expected reductions in carbon dioxide and methane (as equivalent CO₂) emissions to determine the net impact on radiative forcing. In our illustrative future global hydrogen economy scenario (Sect. 2.3), we estimate additional H₂ emissions of between 9 and 95 Tg H₂ yr⁻¹, based on a hydrogen economy supplying 23 % of present-day energy consumption and H₂ leakage rates of 1 % and 10 %. Using a H₂ GWP(100) of 12, this is equivalent to the time-integrated radiative forcing from carbon dioxide emissions of about 110 and 1140 Tg CO₂ yr⁻¹ respectively. Based on the fossil fuel energy sectors that are replaced by hydrogen, a low-carbon method of hydrogen generation and the sector emissions described in Hoesly et al. (2018), we would expect a reduction in carbon dioxide emissions of ∼ 26 000 Tg yr⁻¹. In addition, a further reduction of ∼ 1200 Tg yr⁻¹ equivalent CO₂ emissions would result from expected methane emission reductions (assuming a GWP(100) for methane of 28), resulting in a total reduction of CO₂ equivalent emissions of ∼ 27 200 Tg yr⁻¹. Therefore, in this global scenario the increase in equivalent CO₂ emissions, based on 1 % and 10 % H₂ leakage rates, offsets approximately 0.4 % and 4 % of the total equivalent CO₂ emission reductions respectively. Whilst the benefits from equivalent CO₂ emission reductions significantly outweigh the disbenefits arising from H₂ leakage, they clearly demonstrate the climate importance of controlling H₂ leakage within a hydrogen economy.

A switch to a hydrogen economy would also provide the opportunity to reduce emissions of other gases, which can themselves directly or indirectly affect both climate and air quality. For example, an immediate impact of increased atmospheric H₂ is the reduction of the concentration of OH, the major tropospheric oxidant, which would thus tend to increase the lifetime of methane. Increases in methane would adversely affect climate and also lead to production of tropospheric ozone, impacting both climate and air quality. However, we show that reductions in methane emissions associated with decreased fossil fuel use, along with reductions of CO and NO_x emissions, could lead to a net change in ozone close to (or below) zero with associated climate and air quality benefits, assuming very high H₂ leakage rates (e.g. scenario 2000H2_CH4f_O3pre, Fig. 2j) and a green method of hydrogen production. In the case of no leakage of H₂, emission reductions in methane, CO and NO_x lead to decreases in tropospheric ozone column globally (scenario 500H2, CH4f_O3pre, Fig. 2h). To assess whether the ozone response to changes in emissions associated with a hydrogen economy is influenced by the assumed background climate conditions (i.e. whether the timing of the implementation of a hydrogen economy could be important), we performed a set of additional model simulations using 2045–2055 climatological boundary conditions taken from the CMIP6–SSP2-4.5 scenario. These scenarios assumed a hydrogen economy of the same absolute size and identical leakage rates to our present-day (2014) scenarios. We found no discernible difference between the tropospheric and stratospheric ozone response to changes in H₂ and ozone precursors. Emission reductions, from decreased fossil fuel use, offer similar air-quality benefits in both sets of simulations. To maximise the benefit of a switch to hydrogen, not only should hydrogen leakage be minimised, but the emissions of methane, CO, VOCs and NO_x should also be reduced to the maximum extent possible.