Prather (1994) presented a simple three-equation system to demonstrate that methane perturbations always decay more slowly than methane’s lifetime in the bulk atmosphere. The classic three-equation system considers the following chemical reactions: C H 4 + O H → k 1 m u l t i p l e   s t e p s → C O + p r o d u c t s , R C H 4 = k 1 O H C H 4 C O + O H → k 2   p r o d u c t s ,   R C O = k 2 O H C O O H + X → k 3   p r o d u c t s , R X = k 3 O H X where R i are the reaction rates, k i are the reaction rate constants, and [·] denote concentrations. Only products of interest are named in these equations. Inclusion of the surrogate term X accounts for all additional OH sinks. The time-dependent mass balance equations for each species are then: d C H 4 d t = S C H 4 − R C H 4 d C O d t = S C O + R C H 4 − R C O d O H d t = S O H − R C H 4 − R C O − R X where S i are production rates. To consider hydrogen’s impact on this system, Prather’s three-equation system has been extended to include an additional equation for hydrogen (Bertagni et al., 2022; Warwick et al., 2022; Warwick et al., 2023). The resulting four-equation system considers the following chemical reactions: C H 4 + O H → k 1 m u l t i p l e   s t e p s → C O + p r o d u c t s , R C H 4 = k 1 O H C H 4 H 2 + O H → k 2   p r o d u c t s ,   R H 2 = k 2 O H H 2 C O + O H → k 3   p r o d u c t s , R C O = k 3 O H C O O H + X → k 4   p r o d u c t s , R X = k 4 O H X

The time-dependent mass balance equations for each species are: d C H 4 d t = S C H 4 − R C H 4 − R S , C H 4 d H 2 d t = S H 2 + α R C H 4 − R H 2 − R S , H 2 d C O d t = S C O + R C H 4 − R C O − R S , C O d O H d t = S O H − R C H 4 − R H 2 − R C O − R X

We include soil and stratospheric sinks as first-order loss terms R s , i = k s , i i . Stratospheric and soil sinks for CH₄ and CO are represented by first-order loss rates of 0.015 years⁻¹ and 0.8 years⁻¹ (Szopa et al., 2021). H₂ loss to soils is represented by a loss rate of 0.4 years⁻¹, consistent with the average of current estimates (Warwick et al., 2022; Warwick et al., 2023). The CH₄ and CO soil sinks are both minor relative to their loss to OH, but the H₂ soil sink is the dominant loss term in the hydrogen budget, accounting for about 80% of H₂ loss (Xiao et al., 2007). The hydrogen soil sink is highly uncertain, so we evaluate our results with the current range of H₂ soil sink estimates which correspond with lifetimes against loss to soils of 1.7–2.9 years (Warwick et al., 2022; Warwick et al., 2023). A portion of the CH₄ oxidation produces hydrogen through photolysis of formaldehyde formed as an intermediate, and the term α represents the fractional proportion of CH₄ oxidation that produces H₂ in this way. We set α to 0.59 to correspond with the value obtained for equivalent conditions in our photochemical box model. Reaction of OH with species other than H₂, CH₄, or CO is included in R_x (see Section 2.2). In our photochemical box model, 30% of OH reacts with species included in X, and we set k₄ [X] to 0.25 s⁻¹ so that the 4-equation and full chemistry models have the same partitioning of OH loss between X, CH₄, H₂, and CO. Source term values, steady-state burdens, and rate constants for the four-equation system are shown in Table 1.

Our zero-dimensional photochemical box model uses a chemical mechanism adapted from Daniel and Solomon (1998). Photodissociation reactions are shown in Table 2. Photolysis rates are calculated using the Tropospheric Ultraviolet Visible radiative transfer model (TUV 5.4; Madronich, 1987; Madronich and Weller, 1990) using the two-stream pseudo-spherical option. We calculate the full annual cycle of photodissociation rates at an altitude of 3 km, at a latitude of 15°N, with cloud optical thickness of 9, and with a total ozone column of 250 DU. Our parameters reflect typical conditions in the tropics, the location where the bulk of CH₄ oxidation occurs (Liang et al., 2017).

The chemical reactions considered are shown in Table 3. We calculate rates following Jet Propulsion Laboratory (JPL) 2019 recommendations (Burkholder et al., 2019) at a pressure of 800 mbar and a temperature of 272 K. This temperature was chosen as a suitable representative temperature for OH oxidation of methane following Spivakovsky et al. (2000).

Nonmethane hydrocarbon chemistry is approximated by the chemistry of ethane (Table 3). Direct emission of ethane is set to yield a concentration (∼30 ppb) that leads to a fractional CO production by nonmethane hydrocarbons consistent with model estimates (Holloway et al., 2000). Direct emission of CH₄, CO, NO, and H₂ are set to values which lead to typical clean tropospheric concentrations (Table 4).

The model setup produces an amount of OH that leads to a methane tropospheric photochemical lifetime of 9.7 years and a total methane lifetime of 8.5 years, in general agreement with Szopa et al. (2021). We use the same soil and stratospheric loss rates for CH₄ and CO, and the same soil loss rate for H₂ as in the four-equation system described in Section 2.1. Stratospheric loss of H₂ is negligibly slow (Xiao et al., 2007) and is not included. Wet deposition is simulated by prescribing first-order loss rates of 0.25 days⁻¹ for CH₂O, C₂H₅OOH, CH₃OOH, HNO₃, H₂O₂, and HONO.

We solve the resulting set of differential equations representing photodissociation reactions, chemical reactions, direct emission, and wet and dry deposition using the Newton-Raphson method with fixed 20-min time steps (Brasseur and Solomon, 2005). We discuss sensitivity of our results to temperature, pressure, column ozone, latitude, altitude, NO emission, and strength of the H₂ soil sink in the following section.

Figure 1 shows the difference in methane response between the four-equation system and our full chemical mechanism for a 53 ppb H₂ pulse, about 10% of the typical steady-state H₂ burden. The four-equation system methane perturbation peaks at about 1.9 ppb, and the fuller mechanism methane perturbation peaks at about 1.1 ppb. Using a tropospheric total mass of 4.22 Tg × 10⁹ Tg to convert from mixing ratio to mass (Trenberth and Smith, 2005), the four-equation system peaks at 120 ppt CH₄ per Tg H₂, and our fuller mechanism peaks at 71 ppt CH₄ per Tg H₂. Other studies which performed transient H₂ simulations using comprehensive chemistry schemes found peak methane perturbations of 80 ppt CH₄ per Tg H₂ (Derwent et al., 2020) and 27 ppt per Tg H₂ (Field and Derwent, 2021). The transient behavior of methane due to a hydrogen pulse in our photochemical box model lies between these two previous results, while the four-equation system lies substantially higher.

Relative to our fuller mechanism, the four-equation system overestimates the peak methane perturbation in response to the pulse by about 85%, and the time-integrated methane response by over 100% (123%). The four-equation system overestimates the CH₄ perturbation because it uses a constant value for the source of OH, S_OH, throughout the time integration (see four-equation system defined in Section 2.1). This misses key OH production feedbacks; in effect, CH₄ oxidation over the time of integration affects OH production. For example, in the presence of sufficient NO (e.g., polluted Northern Hemisphere), CH₃ radical oxidation produces more than one HO₂: C H 3 O 2 + N O →   C H 3 O + N O 2   Table  2   R 23 C H 3 O + O 2 →   C H 2 O + H O 2 Table  2   R 29 C H 2 O → h v 2 H O 2 + C O   Table  1   R 10

The resulting HO₂ radicals can cycle back to OH by reacting with NO or O₃, or by self-reacting to form H₂O₂ which can then photodissociate to produce two OH radicals. To a lesser degree, additional ozone produced from the H₂ and CH₄ perturbations also contributes to increased OH production. These OH feedbacks have been documented extensively in previous work (Lu and Khalil, 1993; Daniel and Solomon, 1998; Taraborrelli et al., 2012; Voulgarakis et al., 2013; Holmes, 2018; Lelieveld et al., 2016).

Figure 2 shows the results of sensitivity tests of the following model parameters: temperature (±5 K), strength of the H₂ soil sink (±20%), doubling NO direct emission, increasing total column O₃ to 350 DU, decreasing pressure to 700 h Pa, and increasing the altitude (5 km) and latitude (30°N) at which photodissociation rates are calculated, broadly covering regions where nearly all OH oxidation of methane occurs, as discussed in Liang et al. (2017) and Spivakovsky et al. (2000).

As a check on OH in our model, we calculate the lifetime of methyl chloroform using the rate constant for loss against OH from JPL recommendations (Burkholder et al., 2019) and prescribing a lifetime against stratospheric loss of 40 years so that about 10% of methyl chloroform is lost in the stratosphere (Patra et al., 2021). The resulting lifetime of methyl chloroform in our model is 5.0 years, and the average OH concentration is 9.0 × 10⁵ radicals cm⁻³. From observations, global methyl chloroform lifetimes have been estimated to be 4.8 ± 0.3 years (Prinn et al., 1995), 4.5 ± 0.1 years (Krol et al., 1998), and 5.2 − 0.3 + 0.2 years (Montzka et al., 2000). These lifetimes correspond to inferred OH concentrations of 9.7 ± 0.6 × 10⁵ radicals cm⁻³, 1.07 − 0.17 + 0.09 × 10⁶ radicals cm⁻³, and 1.1 ± 0.2 × 10⁶ radicals cm⁻³, respectively. Within uncertainties, our model’s methyl chloroform lifetime and average OH concentration are broadly consistent with these observationally derived estimates.

Another useful metric to consider is the methane feedback factor (Wild and Prather, 2000). Following the formulation from Fiore et al. (2009), the methane feedback factor f is defined as f = 1 1 − s with the sensitivity coefficient s defined as s = ln τ 1 − ln τ 0   / ln B 1 − ln B 0 where τ 1 is the lifetime of methane in the perturbed run, τ 0 is the bulk lifetime, B 1 is the methane burden from the perturbed run, and B 0 is the burden from the unperturbed run. In the limit of small methane perturbations, this factor describes the proportional increase of steady-state methane burden relative to the proportional increase in methane emission. Effectively, this metric describes how sensitive atmospheric OH is to CH₄ perturbations. We calculate the methane feedback factor for both models by increasing methane emission by 5% for the perturbed runs. Our formulation of the four-equation system has a methane feedback factor of 1.7, while our full chemical mechanism has a methane feedback factor of 1.4 for the baseline model and for the sensitivity tests shown in Figure 2. Previous studies using global models have found methane feedback factors f ≈ 1.2–1.5 (Fiore et al., 2009; Stevenson et al., 2013; Holmes et al., 2013; Holmes, 2018; Thornhill et al., 2021). The Intergovernmental Panel on Climate Change (IPCC) Sixth Assessment Report’s assessed value for f is 1.30 ± 0.07 (Szopa et al., 2021). A higher methane feedback factor corresponds with a longer methane perturbation lifetime, a higher integrated methane response, and therefore a higher H₂ GWP due to the methane perturbation. The four-equation system has a much higher methane feedback factor than our full mechanism because it lacks key chemistry which enhances OH production, as described earlier in this section, thereby dampening the CH₄ response to additions of CH₄ or H₂.

The global warming potential (GWP) of H₂ is the time-integrated radiative forcing indirectly induced by a H₂ pulse relative to the time-integrated radiative forcing due to a pulse of an equivalent mass of reference gas over a chosen time horizon H. The GWP of H₂ relative to CO₂ is defined as: G W P H 2 , C O 2 H = A G W P H 2 H A G W P C O 2 H = ∑ i = C H 4 , O 3 , H 2 O ∫ 0 H E R F i t   d t   ∫ 0 H E R F C O 2 t   d t

The absolute global warming potential (AGWP) of H₂ is the sum of integrated effective radiative forcing (ERF) components due to methane, tropospheric ozone, and stratospheric water vapor.

We additionally calculate the GWP relative to CH₄ as a reference gas: G W P H 2 , C H 4 H = A G W P H 2 H A G W P C H 4 H

We use the IPCC formulation for the AGWP of CO₂: A G W P C O 2 H = A C O 2 ( a 0 H + ∑ i = 1 3 a i τ i 1 − exp   − H τ i where a 0 = 0.2173 , a 1 = 0.224 , a 2 = 0.2824 , a 3 = 0.2763 , τ 1 = 394.4 years, τ 2 = 36.54 years, τ 3 = 4.304 years (Myhre et al., 2013), and A C O 2 = 1.33 ± 0.16 × 10 − 5 W m^-2 ppb^-1 (Forster et al., 2021).

In our model, we directly calculate the AGWP of the CH₄ term from a 53 ppb pulse of H₂ as: A G W P C H 4 H = ∫ 0 H E R F C H 4 t   d t = A C H 4 ∫ 0 H Δ C H 4 d t where A C H 4 is the radiative efficiency of methane (in units of W m⁻² ppb⁻¹) and Δ C H 4 is the methane perturbation due to a pulse of hydrogen. We use the radiative efficiency A_CH4 = 5.7 ± 1.4 × 10⁻⁴ W m⁻² ppb⁻¹ which includes both direct and indirect effects due to ozone and stratospheric water vapor (Forster et al., 2021). We separate the methane-induced GWP by source component using radiative efficiencies of ozone and stratospheric water vapor due to CH₄ of A_CH4,O3 = 1.4 ± 0.7 × 10⁻⁴ W m⁻² ppb⁻¹ and A_CH4,H2O = 0.4 ± 0.4 × 10⁻⁴ W m⁻² ppb⁻¹ (Forster et al., 2021).

To include the GWP components due to ozone and stratospheric water vapor induced by a hydrogen pulse, we use the formulation from Warwick et al. (2022) which parametrizes the ozone and stratospheric water vapor responses due to hydrogen from the UKESM1 model. This formulation is also described in detail and applied by Ocko and Hamburg (2022) in their assessment of the time dependence of hydrogen warming impacts for various emissions scenarios. Relative to CO₂, this amounts to GWP contributions due to ozone and stratospheric water vapor of about 24% and 30% of the total 100-year GWP (with the rest being the methane term discussed above). The time-dependent GWP of H₂ relative to CO₂ is shown in Figure 3A.

To obtain the GWP of H₂ relative to CH₄, the AGWP of H₂ is the same as in the CO₂ calculation, and the AGWP of CH₄ is calculated from a 38 ppb pulse of CH₄, about 2% of the typical steady-state CH₄ burden, using the same A G W P C H 4 H formula as for the H₂ pulse case. We calculate the GWP of H₂ relative to CH₄ from fossil and non-fossil sources. Fossil CH₄ adds fossil CO₂ to the atmosphere, while biogenic CH₄ sources produce CO₂ that was recently removed from the atmosphere. We calculate the GWP of H₂ relative to each of fossil and non-fossil CH₄ by assuming 1 kg of fossil CH₄ generates 2.1 ± 0.7 kg CO₂, and 1 kg of non-fossil CH₄ removes 0.7 ± 0.7 kg CO₂ per kg CH₄ due soil removal and loss of partially oxidized products (Forster et al., 2021).

Figure 3B shows the GWP of H₂ relative to CO₂ and CH₄ from fossil and non-fossil sources. The advantages of considering H₂ GWP relative to CH₄ are that it cancels out the uncertainty due to potential CH₄ lifetime errors as well as errors in its radiative properties, and plateaus after about 15 years. The direct comparison of H₂ to fossil CH₄ should be helpful when considering the tradeoffs of using hydrogen in place of methane for fuel.

Our estimate for the H₂ GWP relative to CO₂ is 28 − 11 + 18 on the 20-year time horizon and 10 − 4 + 7 on the 100-year time horizon. Relative to CH₄, our estimate is 0.36 − 0.07 + 0.15 on any time horizon beyond about 15 years for non-fossil methane. In the case of comparing, for example, the relative climate impact of leakage of hydrogen to the leakage of fossil CH₄, the H₂ GWP relative to CH₄ is 0.35 − 0.06 + 0.13 .

The error in our GWP estimates includes the full range of uncertainty due to the sensitivity cases we presented in Figure 2 and uncertainty due to radiative properties of CH₄, CO₂, H₂O, and O₃ at the 90% confidence level (Forster et al., 2021). For GWPs relative to CH₄, the range includes uncertainty due to atmospheric CO₂ production and loss for fossil and non-fossil CH₄ (Forster et al., 2021). The GWPs broken down by component are shown in Table 5.

Estimates for the CH₄ component of H₂ GWP100 relative to CO₂ in global models range from 2.3 to 6.8 (Table 6). Our estimate for the CH₄ component of GWP100 relative to CO₂ is 4.4, in good agreement with the literature range. Compared with global models, a potential limitation of our photochemical box model is its lack of spatial resolution. Hydrogen’s main loss is through a soil sink that is mostly active in the Northern Hemisphere and displays high latitudinal variance (Paulot et al., 2024). This paper does not address whether hydrogen’s globally averaged GWP is affected by the distribution of loss to soils. Despite these limits, the simple box model is useful to understand how variations in model parameters can affect the calculated hydrogen GWP.

Table 7 shows the CH₄ component of H₂ GWP100 as calculated from the sensitivity cases shown in Figure 2. Box model configurations with shorter methane lifetimes tend to have smaller methane components of GWP100, with the high and low hydrogen soil sink cases as exceptions. Across results from global models, higher methane feedback factors generally correspond with larger methane components of GWP100; our results for the methane component of the GWP lie on the low end among modelling groups shown in Table 6 (see also Supplementary Figure S1).

Using the four-equation system to model different hydrogen fuel implementation scenarios, Bertagni et al. (2022) calculated critical hydrogen emission intensities: the hydrogen leakage rate thresholds above which the atmospheric methane burden would increase despite replacing fossil fuels with hydrogen. If hydrogen leakage rates exceed these critical hydrogen emission intensities, this means there is no CH₄ mitigation from hydrogen fuel use.

Here, we carry out such estimates using the chemical response from our photochemical box model. We follow the same methodology described by Bertagni et al. (2022) to calculate the critical hydrogen emission intensities for 500 Tg per year of green and blue hydrogen production scenarios. The energy from this amount of hydrogen is about 15% of current global fossil fuel energy (Bertagni et al., 2022). In addition to calculating critical hydrogen emission intensities required for methane mitigation, as done by Bertagni et al. (2022), we also calculate critical hydrogen emission intensities for lowering radiative forcing. The first calculation weighs only the atmospheric methane response to the changed hydrogen and methane emissions across scenarios, while the second calculation weighs the GWP tradeoff. Calculating hydrogen leakage rate thresholds weighted by GWP provides a more complete picture of the climate impacts of hydrogen fuel implementation. Critical hydrogen emission intensities are lower when weighing by GWP rather than by change in methane burden because the GWP includes warming impacts from tropospheric ozone and stratospheric water vapor in addition to methane.

Figure 4 shows the critical hydrogen emission intensities for hydrogen scenarios for a range of methane leakage rates during hydrogen production (0% methane emission intensity is green hydrogen, while the rest are blue hydrogen scenarios). For methane mitigation, Bertagni et al. (2022) found critical hydrogen emissions intensities of 9% for green hydrogen, 7% for blue hydrogen with 0.2% methane leakage, 4.5% for blue hydrogen with 0.5% methane leakage, and close to 0% for blue hydrogen with 1% methane leakage. For the same scenarios, we calculate critical hydrogen emissions intensities of 13%, 11%, 7%, and 2%, respectively—i.e., critical hydrogen emissions intensities are larger by about 50% for the more comprehensive chemistry as compared to the four-equation system (see, e.g., the differences in methane response calculated in Figure 1). For these scenarios with respect to global warming potential rather than methane mitigation, the critical hydrogen emissions intensities are 8%, 7%, 5%, and 1%, respectively. The hydrogen emission intensities, weighing either by methane mitigation or global warming potential, cross zero at a methane leakage rate of about 1.2%; blue hydrogen with that leakage rate or greater would not diminish methane burden or global warming.