Edinburgh Research Explorer A multi-model assessment of the Global Warming Potential of hydrogen

With increasing global interest in molecular hydrogen to replace fossil fuels, more attention is being paid to potential leakages of hydrogen into the atmosphere and its environmental consequences. Hydrogen is not directly a greenhouse gas, but its chemical reactions change the abundances of the greenhouse gases methane, ozone, and stratospheric water vapor, as well as aerosols. Here, we use a model ensemble of ﬁ ve global atmospheric chemistry models to estimate the 100-year time-horizon Global Warming Potential (GWP100) of hydrogen. We estimate a hydrogen GWP100 of 11.6 ± 2.8 (one standard deviation). The uncertainty range covers soil uptake, photochemical production of hydrogen, the lifetimes of hydrogen and methane, and the hydroxyl radical feedback on methane and hydrogen. The hydrogen-induced changes are robust across the different models. It will be important to keep hydrogen leakages at a minimum to accomplish the bene ﬁ ts of switching to a hydrogen economy.

W hen hydrogen (H 2 ) is produced, transported, stored, and used, some fraction of the gas will leak to the atmosphere. In the existing value chain, there is very little data on the magnitude of these leakages and how this will evolve in a future growing hydrogen economy. Present-day hydrogen average abundances (mole fraction) are about 530 ppbv as reconstructed by firn air-derived estimates and atmospheric measurements 1,2 . Sources of hydrogen include biomass burning, fossil fuel combustion, biological nitrogen fixation, atmospheric photo-oxidation of methane and volatile organic compounds (VOCs) in the atmosphere, and possibly geological sources 3 . Hydrogen is removed from the atmosphere by biological uptake in soils and atmospheric oxidation by the hydroxyl radical (OH). The largest term, and the largest uncertainty term in the atmospheric hydrogen budget is the soil uptake, which in most studies accounts for 65-85% of the total hydrogen sink [4][5][6] . The atmospheric lifetime of hydrogen, defined as the total atmospheric burden divided by the total sinks, is about 2 years 4 . Any Earth-system perturbation that impacts tropospheric chemistry creates a complex chain of events that alter radiatively active atmospheric species, such as methane, ozone, and aerosols, and hence perturbs Earth's radiative budget. Hydrogen is involved in atmospheric chemical reactions that affect the lifetime and abundances of other gases that have an impact on the climate 7 and is thus such an indirect greenhouse gas.
Four main climate impacts are associated with increased hydrogen levels: (1) a longer methane (CH 4 ) lifetime and hence increased methane abundances, (2) an enhanced production of tropospheric ozone (O 3 ) and changes in stratospheric O 3 , (3) an increased stratospheric water vapor (H 2 O) production, and (4) changes in the production of certain aerosols. The most important reaction driving these impacts is the destruction of hydrogen by OH: OH is the most important and powerful oxidant in the atmosphere 8,9 . Oxidation by OH is the major sink of hydrogen, methane, and other compounds in the atmosphere. The levels of OH in the atmosphere depend on other gases, most notably methane, nitrogen oxides (NOx), carbon monoxide (CO), and VOCs. All these processes are highly coupled and lead to chemical feedback processes 10 . Oxidation of methane and VOCs by OH also provides an important atmospheric source of hydrogen. Similarly, the change in chemistry caused by both the uptake of OH and production of H in reaction (1) can lead to changes in ozone. In the troposphere the production of water vapor is negligible compared to the natural water cycle 11 , but in the stratosphere, this reaction can affect the water vapor levels. Finally, changes in OH, ozone, and other oxidants may also affect the formation of particles, especially sulfate, secondary organic aerosols, and nitrate, and alter their size distributions 12 Global Warming Potentials (GWPs) allow comparisons of the global warming impacts of different gases, and hence hydrogen against other technologies. GWP is a useful metric, facilitating, for example, comparison of different implementations of a future hydrogen economy with other technologies that attempt to mitigate climate change 13 . A review on available literature on global warming from methane and tropospheric ozone from hydrogen was conducted, revealing GWP100 values of 0-9.8 with a central value of 4.3 with 95% confidence (Derwent (2018)) 14 . This means that over a 100-year period, a pulse emission of 1 kg of hydrogen leads to as much global warming as 4.3 kg of CO 2. A follow-up study with the STOCHEM-CRI model found a GWP100 value of 5 ± 1 15 . Recent studies have made a step forward in the research of the climate effect of hydrogen emissions by including the combined effects in the troposphere and stratosphere using global 3D models. Using the UKESM1 model with prescribed hydrogen surface concentrations, a much higher hydrogen GWP100 of 12 ± 6 was estimated 16 . Perturbing hydrogen emissions in a coupled model (GFDL-AM4.1) an indirect radiative forcing of 0.84 mWm −2 per Tg yr −1 hydrogen emitted was estimated; mostly due to longer methane lifetime and increased stratospheric water vapor production 17 . This indirect forcing translated into a GWP100 of 12.8 ± 5.2 (90% confidence interval) 13 . A new analysis with the 2D model TROPOS, which reconciled the values of 2D and 3D models found a GWP100 of 8 ± 2 with 95% confidence 18 .
Here, we present new estimates of the GWP100 of hydrogen, using five global atmospheric chemistry models (GFDL, OsloCTM, INCA, WACCM, and UKCA). We calculate GWP100 using a steady-state perturbation approach that takes all effects of atmospheric hydrogen into account in a comprehensive way. The models are forced by present-day (2010) hydrogen and methane surface concentrations. As idealized experiments, we perturb the models with 10% increase in hydrogen and methane surface concentrations, respectively, and combine the results to derive the full chemical response of perturbing hydrogen emissions. This approach allows us to assess separate uncertainties in the terms used to derive the GWP, and thus give us a better estimate of the overall uncertainty. Two models (OsloCTM and GFDL) also contribute simulations driven by hydrogen emissions and deposition fluxes instead of fixed surface concentrations of hydrogen and give similar results, confirming our understanding of chemical feedbacks that occur in emissions-driven simulations.

Results
The hydrogen budget. Figure 1 shows the hydrogen budget as simulated by six models (GFDL, OsloCTM, INCA, UCI, UKCA, and WACCM). The burden of hydrogen is close to 200 Tg (a) and the total lifetime of hydrogen is 2.4 years (1.9-2.7 years model range) (b). The chemical production and loss of hydrogen in the atmosphere are shown in (c). There is a very good agreement between the models on the chemistry of hydrogen, despite the models being quite different. One exception is WACCM, which has lower formaldehyde levels in the troposphere compared to the other models. In our model set-up for prescribed surface concentrations, the amount of hydrogen in the atmosphere (a) is determined by the surface concentrations, the chemical production, and loss, but does not get affected by the soil sink or emissions. The models diagnose their own soil sink (see Methods), and this is used together with atmospheric loss, to calculate the lifetime in the models. The corresponding emissions can be estimated as the total sink and the total production must balance. In d) the soil sink and estimated emissions are shown. The soil sink values have a large range between 44-73 Tg yr −1 , while the estimated emissions range between 29-68 Tg yr −1 . The total lifetime is calculated by the burden divided by the soil sink (as diagnosed in the models) and atmospheric loss, which is equal to the total loss, as it will be in steady state.
Atmospheric composition changes. To investigate the impact of hydrogen on the atmospheric composition, we have perturbed the models (relative to 2010 levels) with a 10% increase in hydrogen surface concentrations (PERT_H2). Because methane is also fixed at the surface we perform an additional experiment, a 10% increase in methane at the surface (PERT_CH4) to indirectly diagnose the impact of hydrogen via methane changes on the atmospheric composition.
By increasing the hydrogen surface concentration by 10%, the models show consistent changes in atmospheric composition. Figure 2 shows the zonal annual mean changes due to hydrogen in OH, methane, water vapor, and ozone from the surface up to 1 hPa.
Methane-induced changes were included by scaling results from the methane experiment to the expected methane change in the hydrogen experiment (see Supplementary Fig. 2).
In the troposphere, all models show a decrease in OH resulting in increased levels and longer lifetime of methane. The reduction in OH also causes an increase in HO 2 , leading to increased production of tropospheric ozone. Previous studies have shown a large spread in effects on ozone layer depletion, ranging from up to 10% in specific areas and under extreme adaption and leakage scenarios [19][20][21] to minor changes including small increases 22,23 . Here, in the lower stratosphere at levels between 100 and 10 hPa, the models all show only small changes to the ozone layer (Fig. 2).
The oxidation of hydrogen higher up in the dry stratosphere increases water vapor concentrations there due to a much longer lifetime than in the troposphere, causing a cooling effect in the stratosphere and an increased greenhouse effect 24 .
Radiative forcing. Figure 3 illustrates and summarizes the radiative forcing effects of increasing hydrogen. Each bar shows global mean effective radiative forcing (ERF) due to one Tg flux of hydrogen, for each model. Here, the model results are normalized by the change in the hydrogen flux in the perturbed simulations (the emissions needed to sustain the concentration perturbation and hence include chemical feedbacks) and the contribution to ERF via methane changes added (see Method and Supplementary Table 3).
Increased hydrogen in the atmosphere leads to higher methane levels. This gives an ERF from methane (green bars) of 0.46 mW m −2 (0.39-0.51 mW m −2 model range) per flux hydrogen (Tg yr −1 ). The methane increase changes both ozone and stratospheric water vapor in addition to the direct changes in ozone and stratospheric water vapor by increased hydrogen. Both these effects are included in Fig. 3. An increase in methane increases its own lifetime 25,26 . The simulated methane feedback factor ranges from 1.36 to 1.55 (Supplementary Table 2) which is a bit high compared to other multi-model studies 1.34 ± 0.06 26 and 1.30 ± 0.07 27 . Our estimated hydrogen feedback is slightly negative, and the factor ranges from 0.95 to 1.0 (Table S01). These feedbacks are also included in Fig. 3. Fig. 1 The hydrogen budget. The hydrogen budget terms for all the models; a total burden (in Tg), b total lifetime and atmospheric lifetime (in years), c atmospheric chemical loss and atmospheric chemical production (in Tg yr −1 ), and d soil sink and estimated emissions (soil sink + atm. loss -atm. production) (in Tg yr −1 ). The model mean is based on the six models with prescribed surface concentrations of hydrogen (left side). The hydrogen budget is also presented in Supplementary Table 1.
We have calculated the ozone radiative forcing by using a monthly three-dimensional kernel for ozone radiative forcing 28 in the two perturbation experiments relative to the control run. This results in a total ozone ERF of 0.40 (0.28-0.48) mW m -2 per Tg yr −1 of hydrogen flux (yellow bars).
In the stratosphere, hydrogen increases water vapor production. We have calculated the radiative forcing for stratospheric water vapor offline, using separate radiative transfer schemes for longwave and shortwave radiation, including stratospheric temperature adjustments. This results in a stratospheric water vapor ERF of 0.19 (0.08-0.29) mW m −2 per Tg yr −1 of hydrogen flux (purple bars).
Sulfate (SO 4 ) aerosols are formed through gas phase oxidation of SO 2 by OH, but also through aqueous phase reaction of SO 2 with ozone and H 2 O 2 (involving many different chemical species) 29,30 . It is therefore not obvious that a reduced OH level from hydrogen emissions leads to less sulfate aerosols. INCA and OsloCTM have a general reduction in sulfate from hydrogen emission whereas the sign of global sulfate burden varies in the GFDL model. The changes in aerosols give rise to additional radiative forcing via aerosol-radiation-cloud interactions. We have therefore calculated the direct and first indirect radiative forcing of aerosols for the models where this was possible, shown as red bars in Fig. 3. In OsloCTM the global mean burden of sulfate decreases from hydrogen emissions, but aqueous phase reactions lead to more sulfate aerosols at low altitudes leading to a strengthening in the indirect aerosol effect and overall negative forcing. In all model simulations the aerosol forcing is weak (see Fig. 3) and for all experiments, it is weaker in magnitude compared to the change in stratospheric water vapor.
Supplementary Table 3 summarizes the main effects of increasing hydrogen in the atmosphere normalized by 1 Tg flux of hydrogen.
Global warming potential. GWP100 for hydrogen is the ratio of the absolute global warming potential (AGWP100) for hydrogen relative to that for CO 2 . AGWP100 is defined as the timeintegrated radiative forcing to a 1 kg pulse emission of a gas over 100 years. Since hydrogen perturbs many gases, we have summed all forcing contributions, assuming that all perturbations from the initial hydrogen pulse have decayed, so the steady-state calculation will match the integrated transient (see Methods). Figure 4 shows the hydrogen GWP100 values for all models, split into the contribution from methane, ozone, and stratospheric water vapor. The model mean hydrogen GWP100 is 11.6 ± 2.8 (one standard deviation). The largest contribution is from changes in methane (44%), followed by ozone (38%) and stratospheric water vapor (18%). The reported GWP100 does not include any aerosol effects which are consistent with IPCC reporting. The number 11.6 ± 2.8 (one standard deviation) is comparable to Warwick et al. 16 (12 ± 6) and Hauglustaine et al. 13 (12.8 ± 5.2, 90% confidence interval), but twice as high as previously published numbers 14,15 which did not consider changes in the stratosphere. The GWP100 values for each component are provided in Supplementary Table 6. In this study we report GWP100 as it is the official emission metric for United Nations Framework Convention on Climate Change (UNFCCC). We have also calculated other time horizons; GWP20 (37.3 ± 15.1) and GWP500 (3.31 ± 0.98) (Supplementary Table 8). By reducing short-lived climate forcers (SLFCs), such as hydrogen, there is potential to slow down global warming in the next 20-25 years. GWP20 is relevant for shorter time horizons, but it might overestimate the impact of SLCFs over CO 2 , as GWP is an integrated metric, and not an end-point metric 31,32 . Since the main focus of a hydrogen Relative changes (in %) in zonal mean concentrations of OH (first column), CH 4 (second column), H 2 O (third column), and O 3 (fourth column) due to a 10% increase in hydrogen for all models; GFDL (a-d), INCA (e-h), OsloCTM (i-l), UKCA (m-o), WACCM (p-s), GFDL-emi (t-w), and OsloCTM-emi (x-aa). This is done by comparing PERT_H2 (Supplementary Fig. 1) plus the contributions from the change in methane caused by hydrogen in PERT_CH4 (scaled to match the methane change in PERT_H2), to CTRL ( Supplementary Fig. 2). In the GFDL-emi experiment this combination was performed within the experiment directly, which is what is plotted here. The UKCA model used, did not model stratospheric water vapor changes. Black dots mark areas where changes between the control and perturbed run are not significant. Significance is here defined as 95% or more of the points with same latitude and height agree on the sign of the change. Because of different perturbation magnitudes in GFDL-emi, GFDL, and OsloCTMemi, these models have been scaled by 25, 4, and 1.65, respectively, to match the 10% perturbation in the hydrogen surface concentrations in the other models (see Methods). economy is the reduction of CO 2 emissions in order to achieve the Paris Agreement's long-term climate target, the use of a short time-horizon may appear less relevant for hydrogen compared to other SLCFs.
Our model ensemble allows us to extend the assessment of uncertainties. To estimate the uncertainties in the GWP100 value to 2.8, we have performed a standard error propagation using model ensemble values and published literature values (as detailed in Methods and Table 1). The largest uncertainty in the estimate of GWP100 of hydrogen is the soil uptake, which is also the largest term in the hydrogen budget, accounting for as much as 80% of the loss of hydrogen 4,5 . If all other factors were known with full certainty, it alone would contribute an uncertainty of ±2.1 (assuming a standard deviation on the soil sink of ±15 Tg yr −1 based on published soil sink values 5,17 and our model spread). Note, however, that the uncertainty linked to the AGWP of CO 2 is comparable in size-a contribution of ±1.8 to the uncertainty, if it were the sole source of uncertainty. The uncertainty in atmospheric loss of hydrogen due to OH is the third largest term, and accounts for ±0.43 when considered on its own.
Following the same method as for hydrogen, we can calculate GWP100 for methane. Our model ensemble has a methane GWP100 of 26.6 (23-33 model range) (Supplementary Fig. 3), comparable to the IPCC AR6 values GWP100 of 27.0 ± 11 (90% confidence interval) 33 . The model mean methane perturbation lifetime is 10.4 years (Supplementary Table 2) which is lower than the IPCC AR6 assessed value of 11.8 years 33 but within the uncertainty of 1.8 years. A higher perturbation lifetime will enhance the methane GWP value.  The percentage uncertainties are calculated using the same method as Table 7SM.8 of ref. 47 , meaning that each percentage contribution is root square summed. We present 1σ uncertainties (68% confidence interval) in addition to the 90% confidence interval assuming a symmetric soil sink uncertainty. We also present the full absolute uncertainty and the source we used to estimate the uncertainty. Where the model ensemble is listed, we have used the standard deviation of the variable from our model ensemble, otherwise, the literature source for the uncertainty estimate is listed. Similar tables for GWP20 and GWP500 can be found in Supplementary Tables 9 and 10.

Discussion
The largest uncertainty in the estimate of GWP100 of hydrogen is the soil sink. This removal is modulated by soil temperature, soil moisture, and the activity of hydrogen-oxidizing bacteria [34][35][36] .
Published global soil removal numbers range from 55 to 88 Tg yr −1 17,37 , and future work should narrow down this large uncertainty. Current global atmospheric models treat this soil removal in a simplified way and need a more sophisticated treatment of soil sink, but it is limited by the resolution. One way to narrow the uncertainty range in the soil sink would be to get a better estimate of the other sources and sinks of hydrogen. In our model set-up, the models diagnose their own soil sink to be used in estimating the lifetime. A high soil sink will be compensated for by a large estimate of hydrogen emissions, as is apparent in Fig. 1 for WACCM (73 Tg yr −1 ), which is outside the range of bottom-up estimates (53.5-60 Tg yr −1 ) ( Table 1 in ref. 17 Fig. 1c (in which there is a good model agreement), the mean soil sink would be 59 Tg yr −1 . This is close to the model mean soil sink (57 Tg yr −1 ). Since the models employ different soil sink schemes, we have also estimated the GWP100 of hydrogen by using this harmonized soil sink value of 59 Tg yr −1 in all models. This leads to a GWP100 of hydrogen of 11.4 ± 2.8 (Supplementary Tables 7  and 11) close to the model mean of 11.6 ± 2.8 with the models' own soil sinks.
Our results are derived by perturbing the surface concentration of hydrogen. To assess the sensitivity of our results to this setup, we also perform emission perturbations using the GFDL and OsloCTM models. In GFDL-emi, GFDL has been perturbed with 200 Tg yr −1 anthropogenic hydrogen emissions, fully coupled for 50 years (as described in 17 ). We have also run OsloCTM with the same emissions files as used in GFDL (OsloCTM-emi) in the control run, but with a smaller perturbation (14 Tg yr −1 ) in the perturbed run. The emission-driven runs are also a test on the location of the perturbation, as anthropogenic emissions are perturbed with most of the emissions occurring in the NH, compared to the surface hydrogen concentration perturbation, which means that a lower radiative forcing is expected since the lifetime of hydrogen is shorter in the NH. We have tested the GFDL model with various strengths of perturbations, and there is linearity in the perturbations. For both models, atmospheric changes and GWP of hydrogen are very similar between the hydrogen emission-driven and the hydrogen concentrationdriven models. Future work should test for emission location dependency, especially NH vs. SH, aviation altitude, and ocean vs. land.
We have calculated GWP100 by assuming that all perturbations from the initial hydrogen pulse have decayed, so the steadystate calculation will match the integrated transient. A major advantage of assembling the GWP from the component terms is that each term represents a fundamental response of atmospheric chemistry to the hydrogen and methane perturbations, and we can compare these across models. Further, from the model range in individual terms we can assess and identify uncertainties in each, allowing us to propagate an uncertainty in GWP that is more than just the range in model GWPs. A notable disadvantage in calculating successive perturbations in this way is ensuring that the higher-order terms become small, and that we have identified all the major chemical couplings. Thus, it will be important to do atmospheric chemistry-model calculations as full emissionsdriven calculations for methane and to ensure that our separation of the problem into two constrained calculations has not missed any important chemical couplings.
We have run the models with present-day background levels of NOx, methane, and CO. These levels are likely to change in the future and may also influence the climate impacts of hydrogen, making it stronger or weaker. For instance, Warwick et al. 16 finds that a reduction in CO, NOx, and VOC emissions alongside increases in hydrogen leads to smaller increase in methane due to a reduction in OH, but that the total response is complex and strongly scenario dependent. Future work should explore these dependencies further, including changes in co-emitted species and air pollution by switching to a hydrogen economy.
By using an ensemble of models with detailed forcing calculations, taking all effects into account including the stratosphere, we have found that the GWP100 of hydrogen is twice as high as earlier estimates, and comparable to recent estimates 13,16 . Our estimated GWP100 of hydrogen is 39% of that to fossil fuel methane GWP100 in AR6 33 . This emission metric can be used in various mitigation policy decisions, by comparing different GHG reduction measures, or life cycle analysis. The potential benefit of switching to a hydrogen economy will depend on hydrogen leakage rates and potential reductions in CO 2 emissions and coemitted species. Various assumptions and estimates of hydrogen leakages vary from 1-10% depending on how hydrogen is being produced, transported, and stored 22,[38][39][40] . The little data that exists on hydrogen leakages comes from assessments and simulations rather than direct measures. Hydrogen leakages will lead to global warming, and it will be important to keep these leakages at a minimum to accomplish the benefits of switching to a hydrogen economy. Therefore, it will be important to develop instruments for leakage detection for monitoring, so we can get a better picture of how large these leakages are today.

Methods
Model simulations. For calculations of GWP from steady state and to investigate changes in atmospheric composition due to increased hydrogen, we ran three sets of simulations: CTRL: A control simulation with present-day atmospheric composition, where hydrogen and methane surface concentrations (mole fraction) are fixed. PERT_H2: As CTRL but the surface concentration of hydrogen is increased by +10%. PERT_CH4: As CTRL but the surface concentration of methane is increased by +10%.
From the PERT_H2 simulation, the change in atmospheric composition and radiative forcing directly due to the increase in hydrogen can be calculated. As the methane concentration is fixed at the surface in these simulations, we need the additional methane experiment, PERT_CH4, to calculate the change in atmospheric composition and radiative forcing due to the change in the methane lifetime in PERT_H2.
Five models (see model description below and Table 2) have run these simulations with prescribed hydrogen and methane concentrations at the surface. Above this model layer, hydrogen and methane are treated as chemically reactive tracers in the atmosphere. Additional simulations were also performed with two models that were run with hydrogen emissions instead of fixed surface concentration fields. The methane concentrations were fixed at the surface. GFDLemi was perturbed with an extra 200 Tg yr −1 anthropogenic hydrogen emissions and ran for 50 years. The methane surface concentration was increased based on the impact of hydrogen onto OH as estimated from a pure hydrogen perturbation run as described in 17 . OsloCTM-emi ran with the same emissions files as used in GFDL-emi, but with a smaller perturbation of 14 Tg yr −1 (a doubling of the anthropogenic emissions). Note that for OsloCTM-emi, the hydrogen at the surface, and hence also burden, is slightly larger in CTRL compared to the prescribed fields used in the concentration-driven run. The increase in methane Model setup. The CTRL simulation uses present-day hydrogen surface abundances. The hydrogen surface concentration field used by all models is on monthly resolution and on a zonally averaged one-degree resolution. The field was generated using data from https://www.esrl.noaa.gov/ using data from 67 stations. The open data includes monthly averages between 1989 and 2005. The data availability varies from station to station, so to create the averaged field we used average monthly values for each month for each station. For each month those averages were fed into a LOWESS smoothing function with fractional smoothing of 0.5 to create a smooth zonal field. Zonal surface field of methane is from input4mip 41  The models have either run with meteorological fields or with nudged meteorology. In GFDL, a higher perturbation is used due to weaker nudging. See Table 2 for details of the individual model set up and model description below.
Radiative forcing calculations. Because methane is held fixed in PERT_H2, the experiment with 10% increase in methane surface concentrations (PERT_CH4) allows us to estimate the impact of hydrogen on surface concentrations of methane. From CTRL and PERT_H2 we directly diagnose the loss of methane through the CH 4 + OH reaction, and the lifetime of methane due to OH, combined with a lifetime of 240 years for loss in the stratosphere (excluding OH loss) and 160 years for soil sink 43 . The total lifetime of methane in CTRL and PERT_H2 can be calculated. The methane flux (i.e. the emissions needed to sustain the concentrations) is calculated as the burden divided by the total methane lifetime in each simulation. The difference in the methane flux between PERT_H2 and CTRL is the additional methane flux for a 10% increase in the surface concentration of hydrogen. The difference in the flux between PERT_CH4 and CTRL is the flux needed to sustain the 10% increase in methane, and this flux includes the feedback factor. Both the methane and hydrogen feedback factor can be calculated as the perturbation lifetime divided by the total lifetime. The increased concentration of methane in PERT_H2 that would have occurred (at equilibrium) if it had not been held fixed (e.g. 43 ,) is then calculated by combining the change in surface concentrations of methane per methane flux (ppb (CH 4 ) (Tg(CH 4 ) yr −1 ) −1 ) in PERT_CH4 with the calculated methane flux in PERT_H2. The methane radiative forcing is calculated by the increase in surface concentration converted using a concentration-to-forcing factor of 0.448 mW m −2 ppb −1,44 and adjustment term of −14% 33 Table 5 summarize the changes in atmospheric composition and effective radiative forcing due to hydrogen only (PERT_H2-CTRL), and methane only (PERT_CH4-CTRL), respectively. For total ozone, stratospheric water vapor, and aerosols, radiative forcing is calculated offline for the PERT_H2 and PERT_CH4 relative to CTRL.
The ozone radiative forcing is calculated from the changes in the monthly mean three-dimensional ozone fields in the perturbations compared to the control simulation multiplied by a monthly three-dimensional kernel for ozone radiative forcing 28 . As in IPCC AR6, we do not add any adjustments and use the radiative forcing values as ERF.
For stratospheric water vapor, the radiative forcing is calculated offline consistently for all models using separate radiative transfer schemes for longwave and shortwave radiation 45 . Stratospheric temperature adjustment is included in the radiative transfer calculations. There are limited studies for adjustment of methaneinduced changes in stratospheric water vapor, and the IPCC AR6 presented an adjustment less than 1 with low confidence. We choose to treat the RF values here as ERF.
Direct aerosol effect and indirect aerosol effect is calculated offline 46 . For the indirect aerosol effect only the change in effective radius is simulated. Thus, no rapid adjustment to microphysical properties such as liquid water content or cloud fraction is taken into account. For aerosols we report RF and not ERF.
Calculation of fluxes and budgets. For calculations of AGWP we calculated the hydrogen flux needed to sustain the 10% increase in surface concentrations of hydrogen. As for methane, this is calculated as the difference in the flux (lifetime divided by the burden) in PERT_H2 and CTRL. The models calculate the atmospheric production and loss, but soil sink and emissions are not needed in the simulations, as the hydrogen at the surface is prescribed. The models do diagnose the soil sink (see Model descriptions below), and this loss in addition to the atmospheric loss is used to calculate the lifetime in the models and hence the estimated flux of hydrogen in the simulations.
The radiative forcing in the PERT_H2 experiment can then be normalized by the hydrogen flux (W m -2 (Tg yr −1 ) −1 ). The indirect forcing via methane is added by normalizing the forcing by the methane flux and multiply by the methane flux per hydrogen flux (Tg CH 4 (Tg H 2 ) −1 ) from the PERT_H2 experiment.
GWP calculations. GWP H2 is the ratio of the absolute global warming potential (AGWP) for hydrogen relative to that for CO 2 . AGWP is defined as the timeintegrated radiative forcing of a 1 kg pulse emission of that gas over a given time horizon. For use in the United Nations Framework Convention on Climate Change (UNFCCC) the GWP is integrated over time horizons of 20 (GWP20), 100 (GWP100), and 500 (GWP500) years. In this study, we focus on GWP100 (values of GWP20 and GWP500 are found in Supplementary).
The GWP100 of hydrogen (following 8.SM.6 from IPCC AR5 [Myhre et al. 2013]) 47 is defined as: As in AR6, we use the effective radiative forcing (ERF) instead of RF. The AGWP100 of CO 2 is 0.0895 10 −12 W m −2 kg −1 yr from IPCC AR6 Table 7.SM.6 48 . The ERF of hydrogen is the sum of the effective radiative forcing contributions by methane, ozone, and stratospheric water vapor.
We rewrite the AGWP from equation (2) to account for results for a run with a pulse emission of mass m em i given in units of kg yr −1 .
Hydrogen does not have a direct radiative forcing, but it has various indirect forcing effects (i.e., methane, ozone, and stratospheric water vapor). These forcing terms (i) are summed: 49 Where ERF iFC t ð Þ is the effective radiative forcing for compound FC caused by the pulse of species i in the model.
Following the transient decay of emissions in a CTM can be difficult as one must keep accurate track of the rise and fall of the chemical species perturbed by hydrogen. The methane chemical mode has a lifetime of 7-14 years, which means integration needs to be performed over 60 years to accurately account for its rise and decay.
The integral over these transients in all species has been demonstrated to be equal to the steady-state pattern of perturbations from constant emissions scaled by the lifetime of the emitted species 50,51 . This is because when a constant emission perturbation has reached steady state, it contains perturbations from pulses of all ages at one time. The pulse integration result is equal to a simple comparison between forcing changes per burden change in steady state runs with control and constant emission changes scaled with the emitted species lifetime.
Where ΔBurden H 2 is the burden change between a perturbed and a control run that have both reached steady state, and τ H 2 is the lifetime in the control run multiplied by a feedback factor to account for the shift in hydrogen lifetime induced by hydrogen changes. In practice, we use a slightly different, but equivalent expression for the AGWP to calculate it from the experiments. Including the hydrogen feedback factor in the lifetime, means that the AGWP is equal to the steady-state radiative forcing (W m −2 ) divided by flux changes (Tg H 2 yr −1 ) needed to sustain the perturbation.
Where Φ ssÀPerturbed i and Φ ssÀControl i are the steady state fluxes of hydrogen in the PERT_H2 and CTRL runs. In our setup, changes due to methane changes are diagnosed in a separate experiment, which leads to a scaling for methane changes. These correspond to radiative forcing changes between the PERT_H2 and CTRL and scaled PERT_CH4 and CTRL experiments respectively.
For uncertainty calculation, we expand the equation in more detail and rewrite the flux changes as burden changes divided by hydrogen lifetime: Δ denotes changes between the control and perturbed experiments where changes in between PERT_CH4 and CTRL have been scaled to match the expected changes in methane between PERT_H2 and CTRL, if the methane had not been fixed at the surface. r CH 4 is the radiative efficiency of methane including rapid adjustments, ff CH 4 is the methane feedback factor, and C CH 4 is the scaled concentration change of methane.
As mentioned above, the longest-lived, climate-relevant perturbations induced by molecular hydrogen emissions, will be those tied to the methane perturbation with an e-fold time of about 12 yr. Thus, for AGWP100, we can assume that all perturbations from the initial hydrogen pulse will have decayed, and the steadystate calculation will match the integrated transient. For AGWP20, however, perturbations tied to methane will still exist after 20 years, and the steady-state approximation is corrected using the time scale of the methane perturbation (see Supplementary Table 8).
Uncertainty calculations. To estimate uncertainties in the GWP value, we use equation (7) and perform a standard error propagation. We are specifically interested in errors that are due to uncertainties in the soil sink and atmospheric loss, hence we take advantage of the fact that lifetime is burden divided by loss to rewrite: Where we have rewritten the lifetime as the hydrogen burden divided by the loss. Since the lifetime used is the perturbation lifetime, we also need to include the hydrogen feedback factor to keep the equivalence. The first fraction in the expression, , is now the inverse of the relative burden change. We propagate the error terms in the following way: where we have expanded the forcing of methane to include uncertainties in the radiative efficiency r CH 4 and feedback factor ff CH 4 . We find a total uncertainty of 2.8. For the terms that do not measure expanded uncertainty, we use the model mean value. For each of the sigma-terms, we use either the standard deviation of the model spread for the model ensemble, or values from literature. Details on the relative contribution to the total uncertainty for each term, and the source for the spread can be found in Table 1 below. For the soil sink specifically, we report a symmetric uncertainty range of ±15 Tg yr −1 around our model mean of 59 Tg yr −1 . This does not completely cover the interval stated in Ehhalt and Rohrer 5 . If we assume that the interval in Ehhalt and Rohrer 5 of 60 þ30 À20 Tg yr −1 is a 90% confidence interval (this is not explicitly stated in the paper), the range is well covered by our range, except the highest soil sink values stemming from top-down estimates 5,6,37 . These high soil sink values cannot be covered by our model ensemble without assuming some unknown source, since their chemical production is far higher than we are able to produce. We feel therefore that reporting mainly 1 sigma confidence intervals, is most appropriate, and that some additional uncertainty towards lower GWP values is likely when reporting larger uncertainty intervals unless the soil sink can be better constrained.
Model description. Table 2 gives an overview of the models and how they are set up, and below is a detailed description of the individual models.
OsloCTM3. OsloCTM3 52 covers the troposphere and the stratosphere, with comprehensive chemistry in both regions. The model is driven by 3-hourly meteorological forecast data by the Open Integrated Forecast System (OpenIFS, cycle 38 revision 1) at the European Centre for Medium-Range Weather Forecasts (ECMWF) and the horizontal resolution is~2.25 × 2.25°with 60 vertical layers ranging from the surface and up to 0.1 hPa. The model includes 174 components (chemical gases and aerosols) and a complex set of photolytic, bi-molecular, trimolecular, and heterogeneous reactions. The model has aerosol modules for sulfate, nitrate, black carbon, primary organic carbon, secondary organic aerosols, mineral dust, and sea salt (Lund et al. 2018) 53 . Based on estimates of natural and anthropogenic emissions, dry-deposition velocities, wet-deposition for watersoluble species, and solved chemical reactions, changes in the atmospheric composition of the components are calculated. For hydrogen the dry deposition or soil sink scheme used are Sanderson et al. 54 that takes into account soil moisture effect on dry-deposition velocities for different vegetation types, no uptake when snow covers the ground and a reduced rate for cold surfaces according to Price et al. 55 . The dry-deposition values are scaled in the concentration-driven run so that the estimated total emission is~32 Tg yr −1 similar to what was used in Paulot et al. 17 . The same deposition scheme and values were used in the emission-driven runs.
WACCM6. WACCM6 56 is the whole atmosphere version of the Community Atmosphere Model version 6 (CAM6) run with 88 vertical levels and 1.875°× 2.5°h orizontal resolution. It fully represents the meridional overturning circulation of the stratosphere, as well as full tropospheric and stratospheric chemistry, with interactive oxidants and ozone. The deposition velocity is computed within the land model (each land type separately). The configuration was run with the Chemistry AeroClim option which is a full chemistry version of the model run with aerosol climatologies. The water vapor in the model was held fixed with respect to the chemistry such that any change in the water vapor mass mixing ratio due to chemistry was not applied to the tracer. This version was modified to run with the chemical tracers in UKCA uncoupled from the radiation such that the radiation was calculated from fixed 2010 trace gas mass mixing ratios and not informed by the chemistry. This full uncoupling of the chemistry and radiation means the driving meteorology in all simulations is identical, so that differences between simulations are purely due to the applied changes in hydrogen or methane. As the water vapor was fixed, the change in water vapor due to changes in hydrogen and methane was not calculated in these experiments.
The hydrogen deposition in UKCA is based on the scheme developed by Sanderson, 2003. Deposition velocity is linear with soil moisture in all regions except C4 surface types where the deposition has a quadratic-log dependence with soil moisture and shrub surface types where deposition velocity does not depend on soil moisture. Chemical species and reactions specific to the middle atmosphere are also included with chemical species belonging to the chlorine and bromine chemistry 66 . For aerosols, the INCA model simulates the distribution of aerosols with anthropogenic sources such as sulfates, nitrates, black carbon (BC), organic carbon (OC), as well as natural aerosols such as sea salt and dust. The heterogeneous reactions on both natural and anthropogenic tropospheric aerosols are included in the model 60,61,67 . Ammonia and nitrates aerosols are considered as described by Hauglustaine et al. 61 . Based on Hauglustaine and Ehhalt 68 , the Net Primary Productivity (NPP) is used to constrain the seasonal and geographical distribution of H2 drydeposition velocity in LMDz-INCA. The NPP is provided to INCA as monthly mean data prepared with the ORCHIDEE land surface model. GFDL-AM4.1. GFDL-AM4.1 is the atmospheric component of the Earth-System Model 4.1 69,70 with a horizontal resolution of~100 km with 49 vertical levels. The model is run with two different configurations. In GFDL-emi the model is run with repeating hydrogen emissions over 50 years using 2010 conditions as described in Paulot et al. 17 . In the perturbation experiment the hydrogen anthropogenic emissions are increased by 200 Tg yr −1 and surface methane concentrations are increased from 1808 ppbv to 2005 ppbv. In GFDL, the model is run with prescribed hydrogen and methane as the other models, with horizontal winds nudged to 6hrwind speeds from NCEP. A greater perturbation is used (+40%) to reduce noise in the stratosphere, where nudging is weak.

Data availability
Data from the modeling runs described in this paper are available in form of netcdf files with https://doi.org/10.11582/2023.00024 in the NIRD research data archive: https:// archive.sigma2.no/pages/public/datasetDetail.jsf?id=10.11582/2023.00024.

Code availability
Custom code for this paper is available in the form of jupyter notebooks as version v0.1.0 of this code https://github.com/ciceroOslo/Hydrogen_GWP on github, released under the Apache-2.0 license.