Global atmospheric hydrogen chemistry and source-sink budget equilibrium simulation with the EMAC v2.55 model

In this study, we use an earth system model with detailed atmospheric chemistry (EMAC v2.55.2) to undertake simulations of hydrogen (H₂) atmospheric dynamics. Extensive global equilibrium simulations were performed with a horizontal resolution of 1.9°. The results of this simulation are compared with observational data from 56 stations in the National Oceanic and Atmospheric Administration (NOAA) Global Monitoring Laboratory (GML) Carbon Cycle Cooperative Global Air Sampling Network. We introduced H₂ sources and sinks, the latter inclusive of a soil uptake scheme, that accounts for bacterial consumption. The model thus accounts for detailed H₂ and methane (CH₄) flux boundary conditions. Results from the EMAC model are accurate and predict the magnitude, amplitude and interhemispheric seasonality of the annual H₂ cycle at most observational stations. Time series comparison of EMAC and observational data produces Pearson correlation coefficients (r) in excess of 0.9 at eight remote stations located in polar regions and on high mid-latitude islands. A further 23 stations yielded correlation coefficients between 0.7–0.9, predominantly located in remote marine stations across all latitudes and also in polar regions. The quality of model predictions (r<0.5, 9 stations) is reduced in anthropogenically highly polluted stations in east Asia and the Mediterranean region and stations impacted by peat fire emissions in Indonesia, as local and incidental emissions are difficult to capture. Our H₂ budget corroborates bottom-up estimates in the literature in terms of source and sink strengths and overall atmospheric burden. By simulating hydroxyl radicals (OH) in the atmosphere leading to a CH₄ lifetime in agreement with observationally constrained estimates, we show that the EMAC model is a capable tool for undertaking high accuracy simulation of H₂ at global scale. Future research applications could target the impact of potentially significant natural and anthropogenic H₂ sources on air quality and climate, reducing uncertainties in the H₂ soil sink and impacts of H₂ release on the future oxidising capacity of the atmosphere.

1 Introduction

H₂ represents an essential energy vector for 2050 net zero decarbonisation targets to be met. Current demand for H₂ equates to approximately 95 Mt per year with existing uses in the refining industry, as well as the chemical industry for production of ammonia, methanol and other chemicals (Hydrogen Council, 2021; International Energy Agency, 2023). H₂ is also used in the direct reduction of iron along with smaller uses in electronics, glassmaking and metal processing (International Energy Agency, 2023). With increased governmental, financial and policy support, demand for H₂ is forecasted to rise to between 430–690 Mt by 2050 (Hydrogen Council, 2021; International Energy Agency, 2023). Achieving this projected level of hydrogen demand can support clean energy use in: (1) hard to abate sectors such as long-haul trucking, shipping and aviation, (2) sectors that require a clean molecule as a chemical feedstock such as for co-firing of natural gas turbines or industrial processes such as steel manufacturing, (3) sectors that require a source of low carbon heat such as for cement and aluminium production or for buildings (Hydrogen Council, 2021; International Energy Agency, 2023). Use of H₂ offers a lot of potential for securing decarbonisation outcomes, provided clean production pathways are prioritised (Hydrogen Council, 2021; International Energy Agency, 2023), carbon capture and storage technologies (if required) work efficiently and at scale (International Energy Agency, 2020), and leakage rates in the H₂ value chain are minimised with sound engineering design (Esquivel-Elizondo et al., 2023; Fan et al., 2022).

Despite these potential advantages for decarbonisation, H₂ has well-documented climate impacts following its release into the atmosphere which represents an important environmental challenge. In terms of climate impacts, H₂ is an indirect greenhouse gas that leads to increases in radiatively active species by increasing (1) CH₄ lifetime due to H₂ competing for the OH sink (2) tropospheric ozone production due to a chain of reactions initiated by the H atom and (3) stratospheric water vapour that enhances radiative forcing (Derwent et al., 2006; Paulot et al., 2021; Ocko and Hamburg, 2022; Warwick et al., 2022, 2023). Since H₂ release affects the oxidising capacity of the atmosphere, it may also lead to changes in the production of sulphate, nitrate and secondary organic aerosols (Sand et al., 2023). Arising from these modelled results, coupled chemistry-climate modelling has a vital role to play before future H₂ infrastructure is installed to ensure that projected increases in H₂ utilisation do not lead to significant adverse consequences for the earth’s atmosphere, air quality and climate.

Simulation of H₂ atmospheric chemistry impacts has attracted significant research attention both in the past few decades (Hauglustaine and Ehhalt, 2002; Schultz et al., 2003; Tromp et al., 2003; Warwick et al., 2004) and at present (Derwent et al., 2020; Paulot et al., 2021, 2024; Warwick et al., 2023) given the likelihood that demand for H₂ usage will grow and potential environmental impacts still require a solution. Previous attempts at simulating hydrogen mixing ratios with coupled chemistry-climate modelling have met variable levels of success at global scale. In this article, we show that the EMAC model is a highly capable tool for capturing (1) the magnitude, amplitude and seasonality of the annual H₂ cycle and (2) the meridional gradients in H₂ mixing ratios. These findings support the conclusion that the EMAC model consistently represents the interplay between the dominating soil sink (i.e. 75 % of all sink terms) and atmospheric photochemical production (i.e. 63 % of all source terms) which is by far the largest source term for H₂ (Table 2).

Table 1 Chemical budgets and lifetimes for H₂ and CH₄. Uncertainties are calculated as the standard deviation of multi-year annual global means. Note that lifetimes are always calculated with respect to global burden (Prather et al., 2012; SPARC, 2013).

2 Materials and methods

In this work, we employ the EMAC model which couples the 5th generation European Centre Hamburg General Circulation Model (ECHAM5; Roeckner et al., 2003, 2004, 2006) to the Modular Earth Submodel System (MESSy) (Jöckel et al., 2006, 2010). Simulations were performed with T63 spectral resolution which produces a spatial resolution of 1.9° (approximately 180–190 km). Simulations were performed with 90 levels up to 80 km above the earth’s surface, encompassing both the lower and middle atmosphere. Chemical reactions in the atmosphere were modelled with version 1 of the Mainz Isoprene Mechanism (MIM1; Pöschl et al., 2000; Jöckel et al., 2006). The model experiment is representative of the present-day (i.e. the year 2020) and uses meteorology for the years 2006–2023, with the first four years used as spin-up time. Flux boundary conditions were employed for both CH₄ and H₂ to overcome issues with the introduction of artificial sources and sinks arising from using Dirichlet boundary conditions with a prescribed mixing ratio at the lower boundary of the atmosphere. H₂ and CH₄ are chemically coupled and have nearly the same chemical lifetime (Table 1). Both compete for the OH radical as a chemical sink, with OH being by far the dominant sink for atmospheric CH₄ (Saunois et al., 2025; see Sect. 4.1 below). Furthermore, atmospheric oxidation of CH₄ is the largest source for H₂ (Ehhalt and Rohrer, 2009). To adequately simulate such a coupled system, the EMAC model uses flux boundary conditions for sources and sinks of both species. To reach a steady-state for the control simulation, the initial conditions for CH₄ and H₂ were obtained from a 15 years long simulation, covering the period 1990–2005. CH₄ was simulated based on the work of Zimmermann et al. (2020), in which emissions of CH₄ and deposition are represented based on the year 2020. Integration of the equations in the simulation uses a time-step of 450 s, and, due to the relatively long lifetime of H₂, precluding diel variability, instantaneous values are outputted every day.

2.1 Emissions

In this work, the goal is to undertake an equilibrium simulation that reaches steady-state mixing ratios representative of present day atmospheric conditions. Therefore, emissions are based on the year 2020, or the closest year prior to 2020, and are repeated for each year which removes any interannual variability. Due to increasing emissions and its long lifetime CH₄ is not in a steady state. Therefore an equilibrium simulation is not fully representative of the atmospheric state in 2020.

For the long-lived tracer CH₄, the a posteriori emissions and the best combination of the rising-CH₄ scenario of Zimmermann et al. (2020) have been applied. In this work, Zimmermann et al. (2020) show that the EMAC model has been efficient in simulating interactive CH₄ mixing ratios over the last two decades. Therein, the model results compare quite well with NOAA and The Advanced Global Atmospheric Gases Experiment (AGAGE) stations and measurements from CARIBIC (Civil Aircraft for the Regular Investigation of the Atmosphere Based on an Instrument Container) flight observations (Brenninkmeijer et al., 2007). Twelve emission categories are considered here, namely, wetlands other than bogs (SWA), animals (ANI), landfills (LAN), rice paddies (RIC), gas production (GAS), shale gas drilling (SHA), bogs (BOG), coal mining (COA), including minor natural sources from oceans, other anthropogenic sources, volcanoes, oil production and offshore traffic, oil-related emissions (OIL), biomass burning (BIB), termites (TER), and biofuel combustion (BFC). Only emissions from bogs, rice fields, wetlands other than bogs, and biomass burning are subject to seasonal variability. Most of the emissions are based on the emission fields of the Global Atmospheric Methane Synthesis (GAMeS) in which processes with similar isotopic characteristics are aggregated into one group (Houweling et al., 1999). For biomass burning, the GAMeS dataset is replaced by the GFEDv4s (Randerson et al., 2017) and is vertically distributed according to a profile suggested in the EDGAR database (Pozzer et al., 2009). The GFEDv4s biomass burning statistics include agricultural waste burning events. A total amount of 601.1 Tg yr⁻¹ of CH₄ is emitted in the model and detailed emissions for each sector can be found in Tables 1 and 3 of Zimmermann et al. (2020), which also describe in detail the emission optimisation process.

H₂ emissions were taken from the RETRO dataset (Lamarque et al., 2010), which was chosen due to its completeness. As for the other sources, we repeated the emissions based on one single year, namely the year 2000. The RETRO database covers the period 1960–2000, and the last year was taken as representative of 2020 emissions, motivated by the stagnation of H₂ emissions in the past few decades (Paulot et al., 2021). A global value of 14.3 Tg yr⁻¹ for anthropogenic emissions is obtained from the RETRO database, as well as 4.8 Tg yr⁻¹ from soil emissions. Biomass burning emissions were obtained from the GFED (Global Fire Emissions Database) database (Giglio et al., 2013), and accounted for 8.35 Tg yr⁻¹. As the RETRO oceanic emissions are outside the range of emissions suggested by the literature (Paulot et al., 2024), these emissions were upscaled to 3 Tg yr⁻¹ so to be within the suggested range (i.e. between 3–6 Tg yr⁻¹). Both the RETRO and GFED databases provide direct estimates of H₂ emissions without relying upon an assumed H${}_{\mathrm{2}}/\mathrm{CO}$ emissions ratio.

For non-GHGs, different emissions were adopted. Anthropogenic sources of short-lived gases are based on CAMS-GLOB-ANTv4.2 and CAMS-GLOB-AIRv1.1 (Granier et al., 2019), and the emissions are estimated with reduction due to lockdowns during the COVID-19 pandemic (Reifenberg et al., 2022). The reduction in the mixing ratio of the OH radical is below 4 % for most of the atmosphere, with the exception of the uppermost troposphere and tropopause region, which is due to reduced flight activity. The small impact on OH is foremost confined to mid-latitudes of the northern hemisphere. Biomass burning emissions are calculated online on a daily basis and rely on dry matter burned from observations and fire type (Kaiser et al., 2012). The emission factors for different tracers and fire types are taken from Andreae (2019) and Akagi et al. (2011). The simulation uses a climatology of the aerosol wet surface density to calculate heterogeneous reactions. It is based on the CMIP5 (Climate Model Intercomparison Project Phase 5) emissions climatology for the years 1996–2005 low S scenario (Righi et al., 2013). The aerosol distribution for radiative forcing calculation is the Tanre climatology (Jöckel et al., 2006). The biogenic emissions of organic species have been compiled following Guenther et al. (1995) and are prescribed in the model in an offline manner (Kerkweg et al., 2006), with the exception of biogenic isoprene and terpenes, for which the emissions are calculated online (Kerkweg et al., 2006).

2.2 Soil sink implementation

We estimate the soil sink using a two-layer soil model (Yonemura et al., 2000; Ehhalt and Rohrer, 2013a; Paulot et al., 2021). H₂ is assumed to diffuse through a dry top layer of soil with no bacterial activity (layer I), which may be covered by an equally inactive layer of snow. In a second layer below the top layer (layer II), the rate of H₂ removal by high-affinity hydrogen-oxidising bacteria (Paulot et al., 2021) depends on both soil temperature and moisture. The resulting deposition rate is parameterised by:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& v_{\mathrm{d}}=\\ & {\displaystyle \frac{\mathrm{1}}{\left(\mathit{\delta }/D_{\mathrm{soil}}\left({\mathit{\theta }}_{\mathrm{w},\mathrm{I}}\right)+{\mathit{\delta }}_{\mathrm{snow}}/D_{\mathrm{snow}}+\mathrm{1}/\sqrt{\left(D_{\mathrm{soil}}\right({\mathit{\theta }}_{\mathrm{w},\mathrm{II}}\left)Ag\right(T\left)f\right({\mathit{\theta }}_{\mathrm{w},\mathrm{II}}/{\mathit{\theta }}_{\mathrm{p}}\left)\right)}\right)}}.\end{array}\end{array}$$

The first two terms in parenthesis in the denominator of Eq. (1) represent diffusion through the inactive soil layer and the snow layer of thickness δ and δ_snow, respectively. The diffusivity of H₂ in soil is given by Millington and Quirk (1959):

$$\begin{array}{}\text{(2)}& D_{\mathrm{soil}}\left({\mathit{\theta }}_{\mathrm{w}}\right)=\left(\right({\mathit{\theta }}_{\mathrm{p}}-{\mathit{\theta }}_{\mathrm{w}}{)}^{\mathrm{3.1}}/{\mathit{\theta }}_{\mathrm{p}}^{\mathrm{2}})D_{\mathrm{air}},\end{array}$$

which depends on the volumetric soil water fraction θ_w and the volumetric soil pore fraction (i.e. porosity) θ_p. The diffusivity of H₂ in snow is given by:

$$\begin{array}{}\text{(3)}& D_{\mathrm{snow}}=\mathrm{0.64}{D}_{\mathrm{air}},\end{array}$$

while the diffusivity of H₂ in air is given by:

$$\begin{array}{}\text{(4)}& D_{\mathrm{air}}={\displaystyle \frac{\mathrm{0.611}\times \mathrm{1013.25}}{p\cdot \left(\right(T+\mathrm{273.15})/\mathrm{273.15}{)}^{\mathrm{1.75}}}},\end{array}$$

where the diffusivity of H₂ in air depends on the air temperature T in °C and the air pressure p in hPa.

The third term in parenthesis in the denominator of Eq. (1) represents H₂ removal in the lower, active layer. The temperature dependence is given by Ehhalt and Rohrer (2011):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}g\left(T\right)& ={\displaystyle \frac{\mathrm{1}}{(\mathrm{1}+\exp [-(T-\mathrm{3.8})/\mathrm{6.7}\left]\right)}}\\ & +{\displaystyle \frac{\mathrm{1}}{(\mathrm{1}+\exp [(T-\mathrm{62.2})/\mathrm{7.1}\left]\right)}}-\mathrm{1},\end{array}\end{array}$$

where T is the soil temperature in °C.

The soil moisture dependence in terms of the water saturation $S={\mathit{\theta }}_{\mathrm{w}}/{\mathit{\theta }}_{\mathrm{p}}$ for eolian sand is given by Ehhalt and Rohrer (2011):

$$\begin{array}{}\text{(6)}& f_{\mathrm{es}}\left(S\right)=\mathrm{0.00936}{\displaystyle \frac{(S-S_{\mathrm{es}}^{*})(\mathrm{1}-S)}{S^{\mathrm{2}}-\mathrm{0.1715}S+\mathrm{0.03144}}},\end{array}$$

where $S_{\mathrm{es}}^{*}=\mathrm{0.02640}$ is the minimum level of water saturation required for microbial activity. For loess loam the soil moisture dependency is given by Ehhalt and Rohrer (2011):

$$\begin{array}{}\text{(7)}& f_{\mathrm{ll}}\left(S\right)=\mathrm{0.01997}{\displaystyle \frac{(S-S_{\mathrm{ll}}^{*})(\mathrm{0.8508}-S)}{S^{\mathrm{2}}-\mathrm{0.7541}S+\mathrm{0.2806}}},\end{array}$$

where $S_{\mathrm{ll}}^{*}=\mathrm{0.05369}$. For a mixture of eolian sand and loess loam we use the weighted mean given by:

$$\begin{array}{}\text{(8)}& f\left(S\right)={\mathit{\phi }}_{\mathrm{sand}}{f}_{\mathrm{es}}\left(S\right)+(\mathrm{1}-{\mathit{\phi }}_{\mathrm{sand}})f_{\mathrm{ll}}\left(S\right),\end{array}$$

where φ_sand is the sand fraction of the soil. The resolution-dependent constant A represents bacterial activity and is adjusted to yield a global mean deposition velocity of 0.033 cm s⁻¹ over land during 2012 to 2015 (Yashiro et al., 2011). Using the 0.25° grid spacing of the ERA5 input data, we obtain A=10.9.

The thickness of the upper soil layer without hydrogenase (i.e. an enzyme in prokaryotes such as bacteria that consume H₂) activity is parametrised by:

$$\begin{array}{}\text{(9)}& {\mathit{\delta }}_{\mathrm{s}}=\mathrm{0.0057}\left(\right({\mathit{\theta }}_{\mathrm{p}}-{\mathit{\theta }}_{\mathrm{w}})/{\mathit{\theta }}_{\mathrm{w}}{)}^{\mathrm{2.5}},\end{array}$$

in sandy loam and

$$\begin{array}{}\text{(10)}& {\mathit{\delta }}_{\mathrm{l}}=\mathrm{0.109}\left(\right({\mathit{\theta }}_{\mathrm{p}}-{\mathit{\theta }}_{\mathrm{w}})/{\mathit{\theta }}_{\mathrm{w}}{)}^{\mathrm{1.8}},\end{array}$$

in loam. Both δ_s and δ_l are expressed in cm. For a mixture of sandy loam and loam with sand fraction φ_sand we use the weighted mean to calculate the soil layer thickness via:

$$\begin{array}{}\text{(11)}& \mathit{\delta }={\mathit{\phi }}_{\mathrm{sand}}{\mathit{\delta }}_{\mathrm{s}}+(\mathrm{1}-{\mathit{\phi }}_{\mathrm{sand}}){\mathit{\delta }}_{\mathrm{l}}.\end{array}$$

The soil water content in the top, dry layer (i.e. θ_wI) is assumed to be the threshold moisture content below which the bacterial activity vanishes, i.e. there is no H₂ uptake in this layer, and is given by:

$$\begin{array}{}\text{(12)}& {\mathit{\theta }}_{\mathrm{wI}}=S^{*}{\mathit{\theta }}_{\mathrm{p}},\end{array}$$

where $S^{*}=S_{\mathrm{es}}^{*}$ is the threshold moisture content for eolian sand and $S^{*}=S_{\mathrm{ll}}^{*}$ for loess loam.

Accordingly, the remaining water within the top 10 cm of soil is between depth δ and 10 cm, resulting in a soil water content for the second layer (i.e. θ_wII) of:

$$\begin{array}{}\text{(13)}& {\mathit{\theta }}_{\mathrm{wII}}={\displaystyle \frac{\mathrm{10}{\mathit{\theta }}_{\mathrm{w}}-\mathit{\delta }{\mathit{\theta }}_{\mathrm{wI}}}{\mathrm{10}-\mathit{\delta }}}.\end{array}$$

We evaluated Eq. (1) using monthly reanalysis data for soil moisture, soil temperature, air pressure, snow depth and snow density with a 0.25° grid spacing from the ERA5 dataset provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) (Hersbach et al., 2020, 2023). The mean soil moisture and soil temperature for the top 10 cm soil layer was obtained by linearly interpolating the ERA5 soil level data. The volumetric soil water content was then uniformly reduced by 6 % (Paulot et al., 2021). Static soil porosity and sand fraction maps with 0.25° grid spacing were obtained from the Land Data Assimilation System (LDAS) (Rodell et al., 2004; GLDAS, 2024).

The global distribution of the H₂ deposition velocity averaged over the 2012 to 2021 period resembles the maps presented by Paulot et al. (2021), but with more pronounced extrema (Fig. S1 in the Supplement). The zonal mean of the H₂ deposition velocity over land also falls within the range of results reported by Paulot et al. (2021). However, the local minimum around 20° N, due to the Sahara, and the maximum around 10° N, where the transition to more humid regions favours soil uptake, are more distinct in this study than in most of those shown by Paulot et al. (2021) (Fig. S2). Thanks to extensive observational records of atmospheric hydrogen concentrations at stations around the world, it is possible to carry out detailed validation of hydrogen deposition based on the resulting concentrations. This method, presented in the following sections, is more informative than local deposition analysis because it is less susceptible to significant variations in surface conditions and is based on a much larger database.

2.3 Observations

EMAC simulations are compared with observational data from 56 stations (with more than 12 monthly values) that form part of the NOAA GML Carbon Cycle Cooperative Global Air Sampling Network (Petron et al., 2024). Data gaps exist at some stations due to the application of quality control procedures, as well as missing data due to impacts of the COVID-19 pandemic. For comparison with the results of the EMAC equilibrium simulation, the observed monthly values have been detrended by subtracting the trend obtained by a sixth order harmonic regression with a linear trend term, while keeping the mid-year values for 2020 fixed.

3 Results

Figure 1 presents a global map of modelled annual mean H₂ mixing ratios as well as the location of observational stations that we use for model inter-comparison purposes. This global map shows the inter-hemispheric gradient for this molecule whereby H₂ mixing ratios are higher in the southern hemisphere compared to the northern hemisphere. This global map also shows the influence of pollution hotspots in Asia, and the influence of biomass burning emissions in central Africa and peat fire emissions in southeast Asia.

Figure 1 Global map of H₂ mixing ratios and the location of observational stations. Model data is averaged over the years 2010–2023 (inclusive) which is representative of the year 2020 for a steady-state simulation, while observational data uses mid-2020 values from a detrending fit using a sixth order harmonic regression technique.

We show time series data for model comparisons with observational data (without gaps in monthly data) from 20 observational stations in Fig. 2. Further comparisons with observational data from another 36 observational stations (with data gaps) are shown in Figs. B1 and B2. 38 of the 56 stations have marine characteristics (Table A1 and Fig. 1), i.e. are located coastally or on islands. 11 stations can be described as polar (latitude $\ge \left|\mathrm{60}\mathit{°}\right|$), and 14 stations are positioned on mountains (i.e. a single station might have several characteristics). These stations usually measure remote, often marine air masses or free tropospheric concentrations. Just 8 stations Hungary (HUN), Mongolia (UUM) , CIBA (CIB, northern central Spain), Shangdianzi (SDZ, China), Wendover (UTA, Utah), Southern Great Plane (SGP, Oklahoma), Israel (WIS), and Gobabeb (NMB, Namibia) fall in none of these previous categories. Inner continental non-mountain stations are sparse, with none in South America and Australia. The two continental stations in Africa (ASK, and NMB) are located in the desert. There are also no monitoring stations in Central Asia. The three Mediterranean stations (CIB, LMP, and WIS) and the continental Chinese (SDZ) and Mongolian (UUM) stations are located in or close to areas with high air pollution levels (Lelieveld et al., 2002; Silver et al., 2025). The Bukit station (BKT, Indonesia) is also regularly impacted by biomass burning and peat fires (Yokelson et al., 2022).

Across these 56 observational stations, the Pearson correlation coefficient (r) exceeds 0.9 for eight remote stations located either in the Arctic or Antarctic regions and on islands at high mid-latitudes. In such cases, the annual cycle of H₂ is modelled excellently in terms of magnitude, amplitude and seasonality. In contrast, nine stations produce correlation coefficients below 0.5. The three stations in the Mediterranean region CIB (−0.17), LMP (0.05) and WIS (−0.29) and the Chinese and the Mongolian stations SDZ (−0.2) and UUM (0.03) show correlation coefficients close to zero or even negative. Negative r values suggest that the EMAC model does not correctly capture the phasing of the annual H₂ cycle which results in pronounced phase mismatches or anti-correlation in Table  and Figs. 2, B1, and B2. These stations are located in regions known for high levels of air pollution. The two tropical coastal stations BKT (0.14) and the Natal station in Brazil (NAT, 0.14) do show rather low r values compared to other tropical stations. They are impacted by biomass burning (i.e. NAT) and peat fire emissions (i.e. BKT). The comparison at the Hungarian station (HUN, 0.39) shows a considerably lower r value compared to the other two central European stations Hohenpeissenberg (HPB, 0.56) and Ochsenkopf (OXK, 0.64). HPB and OXK are mountain stations less susceptible to local mismatches in H₂ soil deposition, as can be seen from the much better matching of model amplitude and phasing with measurements. The Christmas Island station (CHR, 0.43) is the only very remote tropical island station, which has a considerably low r value. At this station, the observational time series is relatively short with hardly any annual cycle. Another 23 stations have correlation coefficients between 0.7–0.9 which demonstrate very good agreement between model and observational data. A number of these stations are located in the mid-latitudes either in the northern or southern hemisphere. The remaining 16 stations produce correlation coefficients between 0.5–0.7 mostly in either remote tropical or mid-latitude regions. Especially the Antarctic stations show an upward trend for atmospheric H₂. The model with its emissions and soil sink repeating the year 2020, cannot capture this feature. This trend coincides with further increasing atmospheric CH₄ concentrations after 2010 following its hiatus of the previous decade (Lan et al., 2026). Table 2 shows that oxidation of atmospheric CH₄ is the largest source term for H₂ (Ehhalt and Rohrer, 2009). To investigate this further in the future, a model simulation with flux boundary conditions for H₂ and CH₄ in transient mode is needed. Overall, the results are very promising and demonstrate the ability of the EMAC model to predict H₂ mixing ratios accurately in most regions of the earth. To provide a visual overview of these results, Fig. 3 provides a global map of Pearson correlation coefficients for comparison of EMAC and observational data.

Figure 2 Time series comparison of observational and EMAC model data for H₂. Results are presented for 20 stations without any gap in monthly data. The sample size for observational data is denoted by n, while r is the Pearson correlation coefficient. Latitudes are denoted by lat.

Figure 3 Pearson correlation coefficient for the intercomparison between EMAC and observational data for H₂. Model data is compared with detrended observational data for the years 2010–2023 (inclusive) to perform this calculation.

Seiler and Conrad (1987)Warneck (1988)Novelli et al. (1999)Hauglustaine and Ehhalt (2002)Sanderson et al. (2003)Rhee et al. (2006)Price et al. (2007)Xiao et al. (2007)Ehhalt and Rohrer (2009)Pieterse et al. (2011)Yashiro et al. (2011)Paulot et al. (2021)Sand et al. (2023)Paulot et al. (2024)

Table 2 Tabulation of the H₂ budget from this study and from literature estimates. Uncertainties are calculated as the standard deviation of multi-year annual means.

^a Up-scaled from 0.5 to 3 Tg yr⁻¹ to match literature recommendations (Paulot et al., 2021); ^b the dry deposition velocity of H₂ (Paulot et al., 2021) has been reduced by 6 % (from the continental
global mean of 0.035 to 0.033 cms⁻¹) to improve simulated H₂ especially in polar latitudes. ^c Multi-model results are presented as a range. VOC = volatile organic compound. IHD = interhemispheric difference.

We also present a plot of the meridional gradient in H₂ in Fig. 4. Overall, meridional gradients in H₂ are captured very well by the EMAC model, notably for stations located in the southern hemisphere, likely because many represent the background atmosphere, whereas many stations in the northern hemisphere are affected by local influences. The model correctly predicts higher H₂ mixing ratios in the southern hemisphere even though the majority of H₂ sources are present in the northern hemisphere. The predicted interhemispheric gradient in H₂ presented here is correct by virtue of the greater soil sink that is present in the northern hemisphere arising from its larger land area (Ehhalt and Rohrer, 2009). Most of the discrepancies between the observed and predicted H₂ mixing ratios exist for a small number of stations within the northern hemisphere mid-latitudes (between 30–60° N) and in the tropics, presumably influenced by local source variability that is insufficiently resolved by our global model. The coverage of many continental land masses by the observational stations is sparse. For South America, Australia, Africa, Central Asia, Siberia and India there are almost no measurements available. This is a problem, especially in validating the soil sink, which can be considered the most uncertain part of the H₂ budget (Paulot et al., 2021).

Figure 4 Meridional gradients in H₂ for EMAC predictions and observational data where stations had more than 12 monthly values. The solid lines were obtained by locally estimated scatterplot smoothing (LOESS; smoothing parameter $\mathrm{2}/\mathrm{3}$, locally linear). Negative latitudes represent south and positive latitudes represent north of the equator. Data are shown for stations with more than 12 monthly values.

For further results, we refer the reader to Appendix B which provides further graphs and tabulated summaries of model performance.

4 Discussion

4.1 Model comparison with observational data

A key feature of the results (Figs. 2, B1, and B2) is the ability of the EMAC model to realistically predict the magnitude, amplitude and seasonality of the annual H₂ cycle at most stations, in unison with that from CH₄ (Zimmermann et al., 2020), with both compounds being modelled with flux boundary conditions and interactive sinks. Promising results are obtained especially for stations that experience remote air masses, for example, in mostly polar regions, which are particularly sensitive to atmospheric transport and chemistry dynamics. The EMAC model results are also quite promising in a range of mid-latitude stations both in the northern and southern hemisphere. In contrast, there are some regions of the globe (Fig. 3) where results are not as promising. For example, the EMAC model predictions are less accurate in the highly anthropogenically polluted Mediterranean region, near the Amazonian region which is impacted by biomass burning emissions, and southeast Asia which is impacted by peat fire emissions. Due to the coarse spatial resolution of 180–190 km and limited information about local and incidental sources, the variability of mixing ratios in these regions is more challenging to capture. This is especially the case for some coastal stations (e.g. NAT, BKT) where the model is limited due to resolution in accurately representing the mixing of marine and continental air. Also the deviations in China and Mongolia (SDZ, UUM) can be partly attributed to a resolution effect. Both stations are located close to strong horizontal gradients in H₂ mixing ratios. The vertical resolution of the lowermost model layers (i.e. thicknesses of 66, 166 and 319 m from the surface upwards) and the representation of the orography influence the comparison. It is important to consider the measurement height relative to the surface and the geographic prominence of the stations. For example, the modelled amplitude of the annual H₂ cycle can be reduced by up to 40 % between the surface and the next model layer for continental stations due to the importance of the soil sink, whereas the H₂ mixing ratios increase with height driven by the strong atmospheric chemical H₂ production. In addition, an interesting model-measurement discrepancy occurs at the Weizmann Institute of Science (WIS) station near the northern Red Sea, where unaccounted for alkane emissions have been attributed to natural seepage from deep water sources (Bourtsoukidis et al., 2020), possibly accompanied by H₂ emissions. Overall, the EMAC model performs favourably at a global scale for simulating H₂ mixing ratios. Comparison with model output from Yashiro et al. (2011) shows that while the EMAC model produces correlation coefficients in excess of 0.7 for over half of the observational stations, the CHASER chemistry-climate model achieves the same result for only one quarter of all observational stations. The annual mean H₂ mixing ratios are well captured by the EMAC model (Table  and Fig. 4), with the exception of CIB ($r=-\mathrm{0.17}$), UUM (r=0.03), and the coastal station Cape Town (CPT) despite its high Pearson correlation coefficient (r=0.73). As mentioned above, the H₂ observational stations do not represent the inner region of continents well. Especially, the small number (i.e. UTA, SGP, UUM) of remote non-mountain, non-hyperarid stations, placed away from highly anthropogenically influenced regions limits its usage in looking in more detail on the parametrisation of the soil sink. At Wendover (UTA, Utah, arid cold climate, r=0.82) and the Southern Great Planes (SGP, Oklahoma, humid subtropical climate, r=0.77), shown in Fig. 2, the model performs well by realistically representing the magnitude, amplitude and seasonality of the observations, indicating that the soil sink parametrisation of H₂ is realistic in those regions. The poor results at the Mongolian station (UUM, arid cold climate, r=0.03), Fig. B1, might be partly attributed to a resolution effect (see above).

To successfully simulate H₂ mixing ratios in the atmosphere, a model needs to correctly resolve the complex interplay between meteorology and chemistry. In terms of chemistry, having the oxidising capacity of the atmosphere represented correctly is a key consideration (Prather and Zhu, 2024). It is largely controlled by the concentration of OH radicals in the troposphere (Lelieveld et al., 2016) and determines the atmospheric lifetime of numerous species including CH₄. The total CH₄ sink is largely dominated by its reaction with OH (see Saunois et al., 2025 for a review). In this sense, CH₄ lifetime is a measure for the total oxidative capacity of the atmosphere. Observational estimates derived from methyl chloroform (CH₃CCl₃) measurements lead to a total atmospheric CH₄ lifetime of 9.1 ± 0.9 years (Prather et al., 2012). Our estimate of 9.3 ± 0.06 years (Table 1) compares well, and is only marginally higher than the range indicated by Prather et al. (2012). This also holds for the global tropospheric chemical CH₄ sink. The supporting information of Prather et al. (2012) states that for the tropospheric CH₄ lifetime based on the reaction with OH (i.e. 11.2 ± 1.3 years) and for the lifetime of CH₄ based on the reaction with tropospheric chlorine (i.e. 200 ± 100 years) yields a combined tropospheric chemical lifetime of 10.6 ± 1.2 years for CH₄. This value compares quite well with our estimate of 10.4 years (Table 1). Model intercomparisons performed by Nicely et al. (2020) suggests that many chemistry models underestimate CH₄ lifetime due to simulating an atmosphere that is overly enriched in OH radicals. Recently, work by Yang et al. (2025) concurs that several atmospheric chemistry models over-predict OH mixing ratios which has implications for CH₄ and H₂ lifetimes. The multi-model estimate from CMIP6 (Coupled-Model Intercomparison Project 6; Collins et al., 2017) and CCMI (Chemistry Climate Model Initiative; Plummer et al., 2021) models used by the Global Carbon Project community for the bottom-up estimates of CH₄ sources yields a total CH₄ lifetime of 8.2 years with a range of 6.8–9.7 years (Saunois et al., 2025). It also shows a large spread, which propagates into increased uncertainties for the derivation of emission budgets. We believe that the EMAC model is realistically capturing the total oxidising capacity of the atmosphere which helps to facilitate high-accuracy prediction of H₂ and CH₄ dynamics.

4.2 Budget and lifetimes

In the atmosphere, CH₄ and H₂ are tracers strongly connected with similar chemical fates. Table 1 shows that CH₄ and H₂ have nearly identical chemical lifetimes both in the troposphere and atmosphere. Furthermore the total sources of H₂ and CH₄, corrected for molecular masses, are very comparable. The biggest difference between these two compounds stems from hydrogen's much larger soil sink which reduces its tropospheric lifetime by approximately a factor of four compared to CH₄.

In Table 2, we compare our H₂ budget derived from EMAC model output with other estimates from the literature. We find that our H₂ budget agrees favourably with bottom-up literature estimates that rely on a combination of emission datasets and model calculations of turnovers and loss rates, but differs from top-down estimates relying on either inverse modelling (Xiao et al., 2007) or analysis of the ²H (i.e. deuterium) budget (Rhee et al., 2006). Our overall budgeting of sources and sinks agrees very well with bottom-up estimates. In addition, our tropospheric H₂ lifetime is in very good agreement with bottom-up estimates. The tropospheric burden is in the upper range of model estimates. Note, that the upper boundary of the tropospheric range is often not clearly defined in the literature, with different definitions e.g. 100 hPa, World Meteorological Organization (WMO), or a climatological tropopause being used. In this study the WMO tropopause definition is used based on a dynamic tropopause in high latitudes and lapse rate being used at low latitudes. The photochemical production is in between the range for bottom-up and top-down estimates (Paulot et al., 2021). These findings suggest that the EMAC model simulates a realistic atmospheric oxidation capacity which is a critical requirement for predicting H₂ mixing ratios well.

Recent work by the United States Geological Survey (Ellis and Gelman, 2024) has developed a simple, zero-dimensional mass balance model coupled with Monte Carlo uncertainty analysis to explore global potential for geological (or gold) H₂ production in the earth's crust. Median modelled estimates of the subsurface H₂ resource are approximately 5.6 × 10⁶ Mt. Ellis and Gelman (2024) estimate that global geological H₂ resources cause an additional global flux of 24 Tg yr⁻¹ from the subsurface to the atmosphere. This is speculative and would add unaccounted H₂ emissions almost of the strength of the current non-photochemical sources. Current knowledge concerning the budget of atmospheric H₂ does not exclude the existence of a large geological H₂ reservoir, and further emphasises the importance of dry deposition for the global atmospheric H₂ budget.

4.3 Suggestions for future applications

Future research efforts in modelling H₂ atmospheric chemistry could build on the current work in three key ways. Firstly, scenarios could be constructed to explore what role geological H₂ (i.e. gold H₂) holds for future atmospheric chemistry. If economically extractable reserves of gold H₂ are found, future utilisation of H₂ would increase well beyond current projections (Hand, 2023; Truche et al., 2024; Ellis and Gelman, 2024). It would therefore be critical to assess the atmospheric chemistry implications of vastly increased H₂ usage. Secondly, it will be critical to assess what impact H₂ use has on the future oxidising capacity of the atmosphere. Clean H₂ use will be associated with significant reductions in the co-emission of criteria pollutants (Galimova et al., 2022) which will influence the formation of atmospheric oxidants such as ozone and OH radicals that constrain CH₄ and H₂ lifetimes (Archibald et al., 2011; Brasseur et al., 1998; Ganzeveld et al., 2010). Thirdly, the H₂ budget is dominated by the land sink (Tables 1 and 2) and future research efforts could help to constrain the important role played by a number of soil properties (e.g. porosity, soil moisture, temperature, and organic carbon content) on terrestrial H₂ uptake (Ehhalt and Rohrer, 2011, 2013a; Paulot et al., 2021; Smith-Downey et al., 2006). The production of H₂ by enzymes in soil (i.e. hydrogenases) could also be considered in a depth-resolved manner as knowledge of the underlying processes improves (Ehhalt and Rohrer, 2013b). Recent coupling of the JSBACH vegetation model to EMAC by Martin et al. (2024) has developed a potential model tool for undertaking on-line H₂ land sink calculations.

5 Conclusions

In this study, we have successfully extended and used the EMAC model to undertake simulations of H₂ atmospheric dynamics, constrained by flux boundary conditions for both H₂ and CH₄. Comparing the EMAC model output with observational data at 56 stations from the NOAA GML Carbon Cycle Cooperative Global Air Sampling Network generally indicates very good agreement at global scale. Excellent results are achieved at remote observational stations and for stations measuring remote and free tropospheric air, suggesting that atmospheric source, sink and transport processes are accurately represented, while model performance is degraded at stations impacted by nearby pollution sources. Our H₂ budget is also in good agreement with bottom-up estimates in the literature. We find that the EMAC model simulates the CH₄ chemical lifetime in excellent agreement with observational estimates, which suggests the model calculates OH radical mixing ratios in a representative manner. The H₂ soil sink, based on a two-layer soil model (Yonemura et al., 2000; Ehhalt and Rohrer, 2013a; Paulot et al., 2021), in combination with monthly ERA5 reanalysis data for soil related parameters has been successfully used by the EMAC model. We conclude that atmosphere chemistry models with such features, capturing the most dominant terms of the atmospheric H₂ budget, should be able to generally simulate station observations of atmospheric hydrogen. This gives confidence that scenario simulations regarding the future H₂ economy will provide reliable estimates of its atmospheric impact.

Appendix A: List of observational stations

Table A1 Observational Stations from the NOAA GML Carbon Cycle Cooperative Global Air Sampling Network. The station code is complemented with p for polar (latitude $\ge \left|\mathrm{60}\mathit{°}\right|$), ma for marine i.e. $\mathrm{coastal}/\mathrm{island}$, and m for mountain for further characterisation.

Appendix B: Additional results

Figure B1 Time series comparison of observational and EMAC model data for H₂. Latitudes are denoted by lat. [Other stations part 1].

Figure B2 Time series comparison of observational and EMAC model data for H₂. Latitudes are denoted by lat. [Other stations part 2].

Figure B3 Pearson correlation coefficient between EMAC CH₄ mixing ratios and observational data. Model data is compared with detrended observational data for the years 2010–2023 (inclusive) to perform these calculations. For a more extensive comparison see Zimmermann et al. (2020).

Table B1 Comparison of mean model and observational H₂ mixing ratios. Δ= Model − Observed, while r denotes the Pearson correlation coefficient. In the case of the BKT station, the EMAC value of one grid cell to the west of the station is used as it is considered more representative.