# Rift-inversion orogens are potential hot spots for natural H2 generation

Frank Zwaan, Sascha Brune, Anne Glerum, Dylan Vasey, John Naliboff,

Gianreto Manatschal, Eric Gaucher

► To cite this version:

Frank Zwaan, Sascha Brune, Anne Glerum, Dylan Vasey, John Naliboff, et al. . Rift-inversion orogens are potential hot spots for natural H2 generation. Science Advances , 2025, 11 (8), pp.adr3418. 10.1126/sciadv.adr3418 . hal-04998684

HAL Id: hal-04998684

https://hal.science/hal-04998684v1

Submitted on 20 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

GEOLOGY

# Rift-inversion orogens are potential hot spots for natural H₂ generation

Frank Zwaan<fnref n="1" />,²*, Sascha Brune<fnref n="1" />,³, Anne C. Glerum<fnref n="1" />, Dylan A. Vasey⁴, John B. Naliboff⁵, Gianreto Manatschal⁶, Eric C. Gaucher⁷

Naturally occurring hydrogen gas (H₂) represents a potential source of clean energy. A promising mechanism for large-scale natural H₂ generation is serpentinization of exhumed mantle material. We study this serpentinization-related H₂ generation during rifting and subsequent rift-inversion orogen development using numerical geodynamic models. Serpentinization-related H₂ generation is best known from rifted margins and spreading ridges. However, because orogens are colder than rift environments, conditions for serpentinization and natural H₂ generation are considerably better in orogenic settings: We find that yearly H₂ generation capacity from serpentinization in the overriding mantle wedge during rift inversion may be up to 20 times larger than during rifting. Moreover, suitable reservoirs and seals required for economic H₂ accumulations to form are readily available in rift-inversion orogens but are likely absent during bulk serpentinization in rift settings. Together with indications of ongoing natural H₂ generation in the Balkans and Pyrenees, our model results provide a first-order motivation for natural H₂ exploration in rift-inversion orogens.

Copyright © 2025 The
Authors, some rights
reserved; exclusive
licensee American
Association for the
Advancement of
Science. No claim to
original U.S.
Government Works.
Distributed under a
Creative Commons
Attribution License 4.0
(CC BY).

INTRODUCTION

A key challenge for the 21st century is the development of sustain-
able energy sources. Molecular hydrogen gas (henceforth “H₂”)
may be one of our best alternatives to hydrocarbon-based fuels, but
present-day synthetic H₂ production is costly (1, 2). However, H₂ is
also generated by a range of natural (bio)chemical processes in the
lithosphere [e.g., (3, 4)], and this “natural H₂” may represent an
excellent source of sustainable energy that, until recently, has been
mostly overlooked e.g., (5–7).

A prominent mechanism for large-scale natural H₂ generation is
the alteration, or serpentinization, of minerals in (ultra)mafic man-
tle rocks when they react with liquid water (5, 8). This mechanism is
very efficient along rifted margins and mid-oceanic ridges that form
during the divergence stage of the Wilson cycle (9), where mantle
rocks are permanently serpentinized close to the ocean floor (10–
13), but industrial extraction of natural H₂ is likely to be unecon-
nomic in these far offshore and deep-water environments. Onshore
exploration is therefore more promising, where serpentinization-
derived natural H₂ is unequivocally sourced from mantle rocks ex-
posed at the surface [e.g., in Oman (14) and New-Caledonia (15)] or
at relatively shallow (several kilometers) depth [e.g., in Albania (16)
and Kosovo, (17)].

However, because serpentinization and natural H₂ generation is
most efficient at temperatures over 200°C (18, 19), and because H₂ is
quickly consumed by a range of (bio)chemical reactions near the
surface (20, 21), the volumes of natural H₂ available from shallow
mantle rocks most likely remain relatively small. Instead, natural H₂
derived from mantle bodies at greater depth (>10 km) could be the
best target for industrial extraction. Fossil evidence of natural H₂

generation has been reported in the European Alps (22–24), where-
as recent research from the Western Pyrenees points to the presence
of natural H₂ at the surface that is potentially being generated from
mantle rocks at depth (21, 25, 26).

The Alps and Pyrenees are rift-inversion orogens, formed by
large-scale inversion of rift basins during the convergence stage of
the Wilson cycle (Fig. 1) (9, 27). They have in common that ultra-
mafic mantle rocks are exhumed [here defined as being brought
above the original depth of the crust-mantle boundary (Moho) by
tectonic and erosional processes] and emplaced in the overriding
plate [e.g., (28, 29)]. This exhumation enhances the likelihood of
mantle material coming into contact with liquid water (of marine,
meteoric, or metamorphic origin), which is a prerequisite for ser-
pentinization to occur, while remaining at favorable temperatures
for efficient hydration. Yet, despite recent efforts (12, 13), our under-
standing of the dynamics and timing of mantle exhumation and
large-scale serpentinization processes in rift-inversion orogens, cru-
cial to assess the associated natural H₂ potential (1), remain limited,
especially because these processes occur at depth.

In this context, geodynamic modeling provides a vital means of
gaining insights into otherwise hidden subsurface processes, and
various researchers have used numerical modeling methods to
study the evolution of rift-inversion orogens [e.g., (27, 30–32)]. Still,
these studies mostly focused on crustal structures, providing very
few constrains on how rift-inversion orogens may facilitate large-
scale natural H₂ generation during their evolution. Here, we present
generalized numerical geodynamic models that include all key first-
order processes required to describe the large-scale thermomechan-
ical evolution of rift-inversion orogens. These models allow us to
assess where and when mantle rocks may be both at favorable tem-
peratures and in contact with water as they are exhumed below the
thinning continental crust during rifting and in the overriding man-
tle wedge during subsequent inversion. Through this approach, we
can quantify how and when tectonic processes enable large-scale
serpentinization and natural H₂ generation linked to exhumed
mantle material over the course of rift-inversion orogen evolution.
In addition, our modelling approach allows us to conceptually link,

<footnote n="1">¹GFZ Helmholtz Centre for Geosciences, Potsdam, Germany. ²Department of Geo-
sciences, University of Fribourg, Fribourg, Switzerland. ³Institute of Geosciences,
University of Potsdam, Potsdam, Germany. ⁴Department of Earth and Climate Sci-
ences, Tufts University, Medford, MA, USA. ⁵Department of Earth and Environmen-
tal Science, New Mexico Institute of Mining and Technology, Socorro, NM, USA.
⁶University of Strasbourg, CNRS, ENGEES, ITES UMR, 7063 Strasbourg, France.
⁷Lavoisier H₂ Geoconsult, Chamonix, France.</footnote>

<footnote>*Corresponding author. Email: frank.zwaan@gfz.de, frank.zwaan@unil.ch</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Example of a model run and definitions used for model analysis. 3D rendering of model B, with 15 Myr of initial rifting, followed by a 20 Myr of tectonic pause before 20 Myr of inversion. The ASPECT geodynamic model is 2D, whereas the FastScape surface processes model, interacting with the top boundary of the ASPECT model, provided a quasi-3D model result. See Materials and Methods for more details on the overall geodynamic modeling approach. (A and B) General model evolution, with the definition of exhumed mantle [mantle material found above the initial base of the crust (i.e., the original Moho)]. (A) Model state after 5 Myr of rifting. (B) Model state toward the end of the model run (t = 54.5 Myr). (C) Definitions used for serpentinization and natural H₂ generation analysis in rift environments and in the overriding mantle wedge in subsequent rift-inversion orogens. Although serpentinization is also known to occur at higher temperatures and pressures, we adopt the 200° to 350°C temperature range (i.e., the "serpentinization window") to conservatively visualize the H₂ kitchen (in yellow), where efficient serpentinization occurs along large active faults, which are assumed to allow for ample water circulation. The adopted serpentinization window is based on the temperature-defined formula for serpentinization from (12), which we use for the actual estimation of serpentinization and associated H₂ generation. This formula follows a bell curve that peaks within the adopted 200° to 350°C temperature range and is near zero outside of it, illustrating the validity of this temperature range for visualization purposes. See the Materials and Methods section and the supplementary data (37) for more details on the serpentinization and natural H₂ generation analysis, as well as the results of various parameter tests.</figcaption></figure>

in time and space, the possible location of natural H₂ generation with potential migration pathways (faults) and potential reservoirs (sediments). Together, these concepts are required to effectively characterize “hydrogen systems” in nature and to identify potential targets for natural H₂ exploration (21, 33).

RESULTS

Here, we use the geodynamic software ASPECT [e.g., (34, 35)] coupled with the surface processes code FastScape [e.g., (36)] to conduct numerical models of rift-inversion orogens (Fig. 1; for more details, see Materials and Methods). We systematically quantify

the impact of initial rifting duration and the occurrence of a post- rift cooling phase on mantle exhumation and associated large-scale H₂ generation during rift-inversion orogen evolution.

Our key model results are presented in Figs. 2 and 3. Throughout each model's tectonic evolution, we quantify as a first step in two- dimensional (2D) section view the following three attributes: the area of exhumed mantle (i.e., the area of mantle material found above the pre-rift Moho depth of 35 km in our models), the area of actively ex- huming mantle (i.e., exhumed mantle that is moving upward), and the area of exhumed mantle material in the “serpentinization win- dow” (i.e., at temperatures between 200° and 350°C), in which we assume that bulk serpentinization could occur through efficient
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Evolution of key rift-inversion orogen styles. Results shown in (A) and (B) are from model A, in (C) to (F) are from model B, and in (G) to (J) from model C, respectively, with snapshot timing for all models as indicated below model C. We indicate the associated serpentinization as a function of rift duration and post-rift cooling before inversion, where insets in (C) and (D) and (G) and (H) are used to highlight important details. Panel (F) serves to highlight the impact of post-rift cooling in model B. Light green: sub-lithospheric mantle; dark green: lithospheric mantle; white: crust; orange: sediments; turquoise: mantle material in the adopted serpentinization window (between the 200° and 350°C isotherms), transparent gray: fault zones; black: active fault zones in mantle; yellow: areas of highly efficient serpentinization (the H₂ kitchen). See also annotation in (A) and Fig. 1C. Limited serpentinization: up to $5 \times 10^9$ kg(km year)$^{-1}$ (during rifting); extensive serpentinization: up to $1 \times 10^{11}$ kg(km year)$^{-1}$ (during inversion). The results from the serpentinization and H₂ potential analysis are presented in Fig. 3.</figcaption></figure>

mantle hydration, if sufficient water would be available [e.g., (18) (Fig. 1C and Fig. 3A to C)]. The extent of the adopted serpentinization window is based on the temperature-dependent serpentinization formula from (12), where serpentinization efficiency follows a bell curve with its peak (i.e., bulk serpentinization) centered in the 200° to 350°C temperature range (Fig. 1C). Note that our serpentinization parametrization is rather conservative because efficient serpentinization can also occur at higher temperatures and pressures (19).

In a second step, we refine the assessment of serpentinization and natural H₂ generation potential in our modelled rift-inversion orogens by following the methodology from (12) that was previously used for rift settings (Fig. 1C). This methodology involves the tracing

of large active mantle fault zones that are assumed to serve as primary conduits for the water circulation required for serpentinization of mantle rocks (Fig. 1C). We subsequently apply the temperature-dependent bell curve formula (12) to calculate the serpentinization capacity (kilogram of serpentinized mantle per year) along these faults, which is defined as per kilometer perpendicular to the model (i.e., along-strike of the rift or rift-inversion orogen in question) (Fig. 3, D to F). Last, we derive the natural H₂ capacity (natural H₂ generation per year) that is defined as mole per kilometer perpendicular to the model by assuming that 1 kg of serpentinized mantle may generate between 100 and 300 mmol of natural H₂ (Fig. 3, D to F) (11, 23).
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Analysis of serpentinization and natural H₂ generation potential during the evolution of key rift-inversion orogen styles. We show results from models A, B and C, respectively (see also Fig. 2). (A to C) Total amount of exhumed mantle (dark blue), amount of actively exhuming mantle, i.e., upward vertical motion >0 m s⁻¹ (turquoise), and mantle area within the adopted serpentinization window (200° to 350°C) over time (red) for each model. The timing of the model insets is indicated with a black arrow (10, 20 and 30 Myr after the start of inversion, respectively). (D to F) Serpentinization capacity and natural H₂ generation capacity over time. Because of the 2D modeling geometry, all values are provided per along-strike unit length. Note that in these plots we assume that 1 kg of serpentinized mantle generates the maximum of 300 mmol of natural H₂ to draw the black curves, whereas the gray area below these curves indicates the possible variation, assuming that 1 kg of serpentinized mantle may generate between 100 and 300 mmol of natural H₂ (11). R: Maximum serpentinization and natural H₂ capacity during rifting [i.e., ~5 × 10⁹ kg(km year)⁻¹ and ~1.5 × 10⁹ kg(km year)⁻¹, respectively]. I: Maximum serpentinization and natural H₂ capacity during inversion [~1 × 10¹¹ kg(km year)⁻¹ and ~3 × 10¹⁰ mol(km year)⁻¹, respectively]. The model results show that serpentinization and associated H₂ generation is enhanced ~20-fold during inversion, provided that sufficient initial rifting is allowed to take place.</figcaption></figure>

Our approach provides a holistic but simplified first-order assessment of natural H₂ potential in rift-inversion orogens; in the Material and Methods section, we present extensive details on model design and important nuances to our model analysis [e.g., the 2D nature of our models, the focus on serpentinization of exhumed mantle material within the 200° to 350°C temperature window and the use of the temperature-determined serpentinization formula from Liu *et al.* (12) for our analysis, the assumption that the water needed for serpentinization is readily available in active fault zones, and the impact of mantle rock composition and water chemistry on natural H₂ generation]. In the supplementary data (37), an overview of all model results produced for this study is presented, along with results of additional parameter tests demonstrating the robustness of our results regarding serpentinization parametrization, and of our assessment of natural H₂ potential during rift-inversion orogen evolution.

## Impact of limited rifting: Model A

Model A represents the impact of a short initial rifting duration (Figs. 2, A and B and 3, A and D). During the initial 5 million years (Myr) of rifting, a symmetric fault-bounded basin form (Fig. 2A). Once rifting ends, shortening ensues right away (Fig. 2B). The rift basin is inverted within 5 Myr, and a largely symmetric orogen consisting of crustal material with a small core of mantle material develops (Fig. 2B). The short rifting phase allows for a small linear increase in exhumed mantle area (until a maximum of ~1 × 10⁹ m²), but almost all mantle material is buried again during subsequent inversion (Figs. 2, A and B and 3A). Consequently, very little mantle material reaches the adopted serpentinization window, and an absence of large faults means very limited water circulation (Fig. 2, A and B). Hence, the serpentinization capacity remains below 1 × 10⁸ kg year⁻¹ per along-strike kilometer (Fig. 3, A and D).
## Impact of longer rifting and post-rift cooling: Models B and C

Models B and C reveal the impact of longer rift durations and a post-rift cooling phase (Figs. 2, C to J and 3, B and C and E and F). After 15 Myr of rifting in model B, mantle material is exhumed near the surface with only a sparse (1 to 4 km) cover of crustal material and sediments (Fig. 2, C and F). The 25 Myr of rifting in model C leads to break up, with only 1 to 4 km of sediment cover remaining (Fig. 2G). Mantle exhumation over time follows a similar linear increase as in model A but is much more prominent (~3.5 × 10⁹ m² and ~7 × 10⁹ m² at the end of rifting in models B and C, respectively Fig. 3, B and C). In both models, the actively exhuming mantle area deviates from the total exhumed mantle trend after 8 Myr (Fig. 3, B and C). The total area of exhumed mantle remains stable during the subsequent 20 Myr of post-rift cooling, apart from a slow gradual decline due to thermal sag and sediment accumulation (Fig. 3, B and C). However, a large part of exhumed mantle material continues to rise due to isostatic balancing of rift shoulder erosion (Fig. 3, B and C). During inversion, most exhumed mantle is buried again, but the formation of an asymmetric orogen causes a wedge of fresh mantle material to be thrust on top of the down-going plate, near the surface, and well within the adopted serpentinization window (Figs. 2, D and E and H to J and 3, B and C).

Large amounts of mantle material are found in the adopted serpentinization window by the end of rifting in models B and C, which further increases during post-rift cooling (Figs. 2F and 3, B and C). Consequently, serpentinization and natural H₂ generation capacities at the end of rifting in models B and C are much higher than in model A (Fig. 3, D to F). However, the maximum serpentinization capacity of ~5 × 10⁹ kg(km year)⁻¹ and the associated maximum natural H₂ generation capacity of ~1.5 × 10⁹ mol(km year)⁻¹ during rifting are similar in both models, suggesting an upper limit imposed by the vertical extent of the adopted serpentinization window (Fig. 3, E and F). Serpentinization drops to zero during post-rift cooling as no strain localizes (and thus no water is expected to circulate in our analysis; see the Materials and Methods) in the serpentinization window (Figs. 2F and 3, E and F).

The subsequent overthrusting of fresh mantle material within a relatively cool orogen focuses strain along mantle fault zones within the adopted serpentinization window, but on a much larger scale than during previous rifting (Fig. 2, C to J). From the start of inversion, serpentinization and natural H₂ generation capacity in both models B and C are around 1 × 10¹¹ kg(km year)⁻¹ and 3 × 10¹⁰ mol(km year)⁻¹, respectively, about 20 times higher than at the end of rifting (Fig. 3, E and F). This trend is only interrupted by the notable drop in serpentinization capacity at *t* ≈ 50 Myr in model B indicating a temporary strain migration away from the serpentinization window (Fig. 3E). These values [also found in other models with sufficient initial rifting; see (37)] indicate that serpentinization and H₂ generation capacity during inversion can be up to 20 times higher than during rifting but may similarly have an upper limit imposed by the vertical extent of the adopted serpentinization window (Figs. 3, E and F and 4).

## DISCUSSION

### Mantle exhumation patterns and natural H₂ potential

Our numerical models show that rift duration is the key factor enabling exhumation and preservation of mantle material in rift-inversion orogens (Figs. 2 to 4). Sufficient rifting ensures that mantle

material is exhumed to begin with, but rifting also facilitates the formation of asymmetric orogens in which mantle rocks may be preserved or further exhumed [as qualitatively observed in previous numerical modeling studies, e.g., (27, 30–32)].

Given sufficient time, rifting can exhume large mantle volumes (Figs. 1 to 4). However, the high syn-rift thermal gradients limit the size of the serpentinization window and serpentinization capacity during rifting [up to ~5 × 10⁹ kg(km year)⁻¹, similar to values calculated for rifting in other studies (10, 12)]. By contrast, our analysis shows that the maximum serpentinization capacity during inversion can be about 20 times higher [~1 × 10¹¹ kg(km year)⁻¹] due to the much larger serpentinization window in a colder orogenic setting and larger fault zones enabling ample water circulation (Figs. 2 to 4). Still, without active deformation, such as during post-rift cooling, water circulation may be reduced or halted, so that serpentinization is limited or stops altogether (Figs. 2F; 3, E and F; and 4). This need for water circulation along (active) faults for mantle rocks to serpentinize can explain the occurrence of relatively fresh mantle outcrops in orogens [e.g., in the Pyrenees (38, 39)], since our models show that mantle material may exhume without being near these large faults.

The serpentinization capacity in our models is convertible to potential natural H₂ generation capacity, which we estimate to be at most ~1.5 × 10⁹ mol(km year)⁻¹ during advanced rifting and up to 20 times higher [~3 × 10¹⁰ mol(km year)⁻¹] during inversion (Figs. 3 and 4). Furthermore, the resulting rift-inversion orogens generally contain large sediment volumes that can act as reservoirs and seals when natural H₂ is generated in bulk during inversion (21). For example, reservoirs may consist of permeable coarse-grained deposits (sands and gravel) formed in proximal environments throughout both rifting and inversion, whereas seals may consist of impermeable fine-grained deposits (clay) that originally formed in distal rift settings before being incorporated in the orogen. In contrast, reservoir-quality sediments are most likely absent in the distal rift environments dominated by fine-grained deposits at the time of bulk H₂ generation during advanced rifting (Figs. 1 to 4).

These reservoirs and seals are crucial for capturing natural H₂ migrating from the serpentinizing mantle (24) and would enable the establishment of a fully fledged orogenic “hydrogen system,” a concept analogous to “petroleum systems” (33). Such a hydrogen system involves an exhumed mantle source rock in the serpentinization window generating natural H₂ (the “H₂ kitchen”) and its subsequent migration to reservoirs with a structure and seal, where the natural H₂ can accumulate for drilling and extraction (Fig. 5) (21). The orogenic hydrogen system concept we propose here is conservative, as it focuses on bulk natural H₂ generation by a H₂ kitchen situated in the exhumed mantle wedge, which we believe provides the best chance for accumulation of economical volumes to form exploitable H₂ fields. However, additional potential, not considered in this study, may lie in contributions from efficient serpentinization and natural H₂ generation at higher temperatures and pressures (19) or natural H₂ generated by less efficient serpentinization elsewhere in the system that may accumulate in considerable volumes over geological time.

Whatever the source of the natural H₂, its preservation during migration and after arrival in a reservoir is crucial for any meaningful accumulation to form. Hence, natural H₂ fields should ideally have reservoir temperatures above the “biotic fringe” (the ~122°C isotherm that indicates the limit of life) above which H₂ is not consumed by microbial activity (20, 40) (Fig. 5). Alternatively, when no
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Schematic summary of the evolution of natural H₂ potential during rifting and subsequent inversion, depending on rift duration and post-rift cooling.</figcaption></figure>

Natural H₂ capacity (H₂ generation per time unit, as indicated by the yellow bubbles in the sections, as well as in the central graph) occurs in the more advanced rifting stages but can be up to 20 times higher during inversion. The assumption that large-scale faulting and associated water circulation is required for serpentinization to occur is highlighted by the absence of natural H₂ generation during post-rift cooling in our analysis. We also specify the likelihood of reservoir (e.g., sands) and seal (e.g., clays) being available to allow for natural H₂ accumulation during the various stages of rift-inversion orogen evolution. Light green: Sub-lithospheric mantle; dark green: lithospheric mantle; white: crust; orange: sediments; turquoise: mantle material in the adopted serpentinization window (between the 200° and 350°C isotherms); transparent gray: fault zones; black: active fault zones in mantle; yellow: areas of highly efficient serpentinization (the H₂ kitchen). See also Fig. 1C for legend. RIO, rift-inversion orogen.

H₂ fields are available, the serpentinizing mantle rocks in the H₂ kitchen themselves could be directly targeted and stimulated (41) for natural H₂ extraction (Fig. 5). This natural H₂ extraction could be combined with other resource exploitation efforts as well (41) (see Materials and Methods).

The key to successful natural H₂ exploration will be the application of a “natural hydrogen system analysis” approach that assesses the spatial and temporal aspects of a natural H₂ system, similar to the “petroleum system analysis” approach commonly used for hydrocarbon exploration (6, 21, 33, 42).

## Promising natural H₂ exploration sites in the Alpine-Himalayan orogenic belt

Our modeling results enable a solid, first-order evaluation of the natural H₂ potential associated with exhumed mantle as part of the overriding plate in rift-inversion orogens. A key result of our study is that relatively cold rift-inversion orogens, such as those in the western Alpine-Himalayan orogenic belt (Fig. 6), provide much better environments for natural H₂ exploration due to their much larger serpentinization window than is the case in relatively hot rift, rifted margin, and mid-oceanic ridge settings.
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Conceptual depiction of a natural H₂ system linked to serpentinization of exhumed mantle in the overriding plate of a rift-inversion orogen. Figure based on the geology of panel (2') from Fig. 4. Water, originating from either the deep or the surface, can flow along (active) fault zones toward those parts of the mantle that are in the (*conservatively defined, see main text for details) serpentinization window (200° to 350°C), where natural H₂ generation is considered most efficient (the H₂ kitchen). The H₂ gas can subsequently migrate toward the surface along fault zones and may accumulate in reservoirs capped by seals of impermeable layers within sedimentary units. Ideally, this natural H₂ accumulation occurs below the “biotic fringe” (the 122°C isotherm that indicates the limit of life) above which H₂ is not consumed by microbial activity (20, 40), so that a conventional natural H₂ field may develop that can be exploited. Alternatively, it would be possible to directly drill into the natural H₂ kitchen and stimulate the mantle source rock to enhance natural H₂ generation (41).</figcaption></figure>

The Pyrenees, which are characterized by a gradient in initial rift basin maturity from limited rifting in the east to hyperextension in the west, represent a potential target area for natural H₂ exploration (29, 43) (Fig. 6, C and D). While no evidence for exhumed mantle exists in the eastern Pyrenees (Fig. 6C) (29, 43), its presence underneath the Mauléon Basin to the west is well established (Fig. 6D) (44–46). Recent studies in this inverted basin showed that natural H₂ generation may be on-going here, compatible with the fact that the well-identified exhumed mantle body is in the serpentinization window (21, 25, 26) (Fig. 6D). Given the conservative maximum natural H₂ generation capacity value of ~3 × 10¹⁰ mol(km year)⁻¹, we obtain for rift-inversion orogens, and assuming a moderate 10 km extent of the exhumed mantle body along-strike, up to 3 × 10¹¹ mol (i.e., 6 × 10⁸ kg) of natural H₂ could be generated below the Mauléon Basin each year. Adopting an average yearly primary energy use of 3.6 × 10⁴ kWh per person in France (47), the caloric value of such a quantity of natural H₂ derived from this mantle body alone (~2 × 10¹¹ kWh) would cover the yearly energy needs of the nearby city of Toulouse and its population of ~500,000 inhabitants. Together with the presence of thick sedimentary deposits, including extensive evaporite layers that form ideal seals, all ingredients may be in place for a functional hydrogen system where natural H₂ can accumulate to volumes considerably larger than the yearly generation rate (Figs. 5 and 6D) (21).

The wider western Alpine-Himalayan orogenic belt contains various rift-inversion orogens that formed after the closure of a number

of small oceanic basins [e.g., (28, 48, 49) (Fig. 6, A, E and F)]. For instance, the Southern Alps contain exhumed mantle material (Ivrea mantle wedge) situated at suitable depths in the retro-wedge (Fig. 6E) (28), which, in combination with indications of fossil natural H₂ generation found in exposed mantle rocks (22–24), makes for a promising natural H₂ exploration target area. Also, the Betics (Ronda peridotite massif) and the Balkans (remnants of the Meliata-Vardar ocean) contain various exhumed mantle bodies that are linked to abiotic gas generation and warrant further exploration (Fig. 6, A and F) (16, 17, 48, 50, 51). On a larger scale, the branches of the Alpine-Himalayan orogenic belt stretching farther eastward into Asia, where various mantle bodies that are exhumed could be of great interest for natural H₂ exploration as well (Fig. 6A) (52).

## Implications and future perspectives

New geo-resources are urgently needed to facilitate the energy transition. To address this need, it is crucial to develop novel concepts and exploration strategies rooted in research-driven approaches. One such concept is the extraction of natural H₂ generated by serpentinization of exhumed mantle material in rift-inversion orogens. Its potential is highlighted by both existing field data [e.g., (16, 17)], and our new modelling results. The current state of natural H₂ exploration may be analogous to that of petroleum exploration before the 1859 Drake oil discovery in Pennsylvania, which triggered the rise of the oil industry (6). Future work is needed to develop robust strategies for natural H₂ exploration in rift-inversion orogens, especially
Tectonic map of western Alpine-Himalayan orogenic belt

<figure><img src="image_8.png" /></figure>

RIOs with limited initial rifting

<figure><img src="image_9.png" /></figure>

Rift-inversion orogens with advanced initial rifting (hyperextension or oceanization) in the western Mediterranean

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Examples of rift-inversion orogens (RIOs) in the western part of the Alpine-Himalayan orogenic belt. (A) Tectonic map of the western Alpine-Himalayan orogenic belt, showing the locations of the sections in (B) to (F), each of which contains a first-order estimation of the in-section area of exhumed mantle material found in the (*conservatively defined, see text for details) serpentinization window (200° to 350°C), assuming a standard geotherm of 25°C km⁻¹ (75). Modified and simplified after Woudloper (76) (used under a CC BY-SA 1.0 license; https://creativecommons.org/licenses/by-sa/1.0/) and Faccenna *et al.* (77). [(B) and (C)] High Atlas and Eastern Pyrenees: Rift-inversion orogen with little natural H₂ potential, as far as serpentinization of mantel material is concerned because no mantle is found in or near the serpentinization window. Modified and simplified after Muñoz (78), with permission of the Geological Society, London, and Arboleya *et al.* (79). [(D) to (F)] Western Pyrenees, Central Alps, and Betics: Rift-inversion orogens with advanced initial rifting that contain exhumed mantle material within the serpentinization window (i.e., the Mauléon body, Ivrea mantle wedge, and Ronda peridotite massif, respectively). Modified and simplified after Mazzoli *et al.* (50), Rosenberg and Kissling (80) (compilation based on sources cited therein), and Wang *et al.* (44), used with permission of the Geological Society, London, and of the Geological Society of America.</figcaption></figure>

given the additional economic opportunities linked to natural H₂ exploitation, such as production of associated CH₄, geothermal energy generation, CO₂ storage, and mineral extraction (23, 41).

These exploration strategies may be informed by our modelling results, which allow the quantification of the mechanisms and timing of tectonic processes that facilitate large-scale serpentinization and natural H₂ generation during rift-inversion orogen evolution (Figs. 2 to 4). Moreover, our modeling work provides first-order insights into the dynamic spatial and temporal links between the location of natural H₂ generation (source/H₂ kitchen), potential natural H₂ migration pathways (faults), and reservoirs (sediments), where natural H₂ may accumulate in orogenic hydrogen systems (Fig. 5). Similar to exploration efforts in petroleum systems, understanding the dynamic evolution of these links will be critical for accurately characterizing hydrogen systems and identifying viable targets for

natural H₂ exploitation as well as associated economic resource extraction opportunities.

## MATERIALS AND METHODS

### Numerical method

We base our modelling approach on previous work by Neuharth *et al.* (36) using the thermo-mechanical geodynamic code ASPECT (Advanced Solver for Planetary Evolution, Convection, and Tectonics) (34, 35, 53–57) coupled with FastScape (36, 58–60) for the inclusion of surface processes (Figs. 1 and 7). We apply a modified version of ASPECT 2.4.0-pre for the coupling with FastScape. ASPECT and FastScape source code and installation details, the custom ASPECT plugins we use, the parameter files used for each model run, the log files of each model run, and the ParaView state files used for model
<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Reference model setup. (A) 3D visualization of the 2D geodynamic model (ASPECT) with the quasi-3D surface processes model (FastScape) on top. Black lines indicate temperature contours, and initial strain (noise) at t = 0 Myr is present in the uppermost 50 km. (I) and (II) indicate the outer and central locations of the yield strength profiles depicted in (B), which have different crustal thicknesses as the upper crust in the model center is 25 km thick, instead of the 20 km at the model edge. (B) Yield strength and temperature profiles along (I) and (II) shown in (A). (C) Plastic (brittle) and viscous weakening intervals applied in our models. Modified after Neuharth et al. (36). UC, upper crust.</figcaption></figure>

analysis are provided in a publicly accessible Zenodo repository (61) (https://doi.org/10.5281/zenodo.14672886).

## Geodynamic modeling (ASPECT)

We apply ASPECT to solve the conversation equation under the extended Boussinesq approximation.

Conservation of momentum is described by the following equation

$$-\nabla \cdot 2\eta\dot{\epsilon}(\mathbf{u}) + \nabla P = \rho g \quad (1)$$

where η is viscosity, ˙ε is the deviatoric strain rate, u is velocity, P is pressure, ρ is density, and g is gravitational acceleration.

Conservation of mass is described as follows

$$\nabla \cdot \mathbf{u} = 0 \quad (2)$$

Conservation of energy is defined as

$$\bar{\rho}C_p \left( \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T \right) - \nabla \cdot k \nabla T = \bar{\rho}H + (2\eta\dot{\epsilon}):\dot{\epsilon} + \alpha \rho T (\mathbf{u} \cdot \mathbf{g}) \quad (3)$$

where ρ̄ is the reference adiabatic density, C_p the specific heat capacity, T the temperature, k the thermal conductivity, H the radiogenic heating, and α the thermal expansivity. The right-hand terms represent radioactive heating, viscous shear heating, and adiabatic heating, respectively. See Table 1 for the used parameter values.

Each compositional field c_i (e.g., representing lithologies such as the upper crust or strain fields) is advected with the calculated velocity field

$$\frac{\partial c_i}{\partial t} + \mathbf{u} \cdot \nabla c_i = q_i \quad (4)$$

where reaction rate q_i is nonzero for the plastic (brittle) and viscous strain fields.

## Landscape evolution modeling (FastScape)

FastScape modifies ASPECT's model domain's surface as a function of stream-power law fluvial erosion, hillslope and marine diffusion,
<table><caption>Table 1. ASPECT parameters.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Units</th><th>Sediment</th><th>Upper crust</th><th>Lower crust</th><th>Lithospheric mantle</th><th>Sub-lithospheric mantle</th></tr></thead><tbody><tr><td>Reference surface density*</td><td>ρ<sub>0</sub></td><td>kg m<sup>-3</sup></td><td>2520</td><td>2700</td><td>2850</td><td>3280</td><td>3300</td></tr><tr><td>Adiabatic surface temperature</td><td>T<sub>AS</sub></td><td>°C</td><td>1284</td><td>1284</td><td>1284</td><td>1284</td><td>1284</td></tr><tr><td>LAB isotherm temperature</td><td>T<sub>LAB</sub></td><td>°C</td><td>1284</td><td>1284</td><td>1284</td><td>1284</td><td>1284</td></tr><tr><td>Thermal expansivity</td><td>α</td><td>K<sup>-1</sup></td><td>3.7 × 10<sup>-5</sup></td><td>2.7 × 10<sup>-5</sup></td><td>2.7 × 10<sup>-5</sup></td><td>3.0 × 10<sup>-5</sup></td><td>3.0 × 10<sup>-5</sup></td></tr><tr><td>Thermal diffusivity</td><td>κ</td><td>m<sup>2</sup> s<sup>-1</sup></td><td>7.28 × 10<sup>-7</sup></td><td>7.72 × 10<sup>-7</sup></td><td>7.31 × 10<sup>-7</sup></td><td>8.38 × 10<sup>-7</sup></td><td>8.33 × 10<sup>-7</sup></td></tr><tr><td>Heat capacity</td><td>C<sub>p</sub></td><td>J kg<sup>-1</sup> K<sup>-1</sup></td><td>1200</td><td>1200</td><td>1200</td><td>1200</td><td>1200</td></tr><tr><td>Heat production</td><td>H</td><td>W m<sup>-3</sup></td><td>1.2 × 10<sup>-6</sup></td><td>1.0 × 10<sup>-6</sup></td><td>0.1 × 10<sup>-6</sup></td><td>0</td><td>0</td></tr><tr><td>Cohesion</td><td>C</td><td>Pa</td><td>5 × 10<sup>6</sup></td><td>5 × 10<sup>6</sup></td><td>5 × 10<sup>6</sup></td><td>5 × 10<sup>6</sup></td><td>5 × 10<sup>6</sup></td></tr><tr><td>Internal friction angle (unweakened)</td><td>φ</td><td>°</td><td>26.56</td><td>26.56</td><td>26.56</td><td>26.56</td><td>26.56</td></tr><tr><td>Plastic strain weakening interval</td><td>–</td><td>–</td><td>[0,1]</td><td>[0,1]</td><td>[0,1]</td><td>[0,1]</td><td>[0,1]</td></tr><tr><td>Plastic strain weakening factor</td><td>Φ<sub>wf</sub></td><td>–</td><td>0.25</td><td>0.25</td><td>0.25</td><td>0.25</td><td>0.25</td></tr><tr><td>Viscous strain weakening interval</td><td>–</td><td>–</td><td>[0,1]</td><td>[0,1]</td><td>[0,1]</td><td>[0,1]</td><td>[0,1]</td></tr><tr><td>Viscous strain weakening factor</td><td>–</td><td>–</td><td>0.25</td><td>0.25</td><td>0.25</td><td>0.25</td><td>1.0</td></tr><tr><td>Creep properties<sup>†</sup></td><td></td><td></td><td>Wet quartzite</td><td>Wet quartzite</td><td>Wet anorthite</td><td>Dry olivine</td><td>Wet olivine</td></tr><tr><td>Stress exponent (dis)</td><td>n</td><td>–</td><td>4.0</td><td>4.0</td><td>3.0</td><td>3.5</td><td>3.5</td></tr><tr><td>Constant prefactor (dis)</td><td>A<sub>dis</sub></td><td>Pa<sup>n</sup>s<sup>-1</sup></td><td>8.57 × 10<sup>-28</sup></td><td>8.57 × 10<sup>-28</sup></td><td>7.13 × 10<sup>-18</sup></td><td>6.52 × 10<sup>-16</sup></td><td>2.12 × 10<sup>-15</sup></td></tr><tr><td>Activation energy (dis)</td><td>E<sub>dis</sub></td><td>J mol<sup>-1</sup></td><td>223 × 10<sup>3</sup></td><td>223 × 10<sup>3</sup></td><td>345 × 10<sup>3</sup></td><td>530 × 10<sup>3</sup></td><td>480 × 10<sup>3</sup></td></tr><tr><td>Activation volume (dis)</td><td>V<sub>dis</sub></td><td>m<sup>3</sup> mol<sup>-1</sup></td><td>0</td><td>0</td><td>38 × 10<sup>-6</sup></td><td>18 × 10<sup>-6</sup></td><td>11 × 10<sup>-6</sup></td></tr><tr><td>Constant prefactor (diff)</td><td>A<sub>diff</sub></td><td>Pa<sup>-1</sup> s<sup>-1</sup></td><td>5.79 × 10<sup>-19</sup></td><td>5.79 × 10<sup>-19</sup></td><td>2.99 × 10<sup>-25</sup></td><td>2.25 × 10<sup>-9</sup></td><td>1.5 × 10<sup>-9</sup></td></tr><tr><td>Activation energy (diff)</td><td>E<sub>diff</sub></td><td>J mol<sup>-1</sup></td><td>223 × 10<sup>3</sup></td><td>223 × 10<sup>3</sup></td><td>159 × 10<sup>3</sup></td><td>375 × 10<sup>3</sup></td><td>335 × 10<sup>3</sup></td></tr><tr><td>Activation volume (diff)</td><td>V<sub>diff</sub></td><td>m<sup>3</sup> mol<sup>-1</sup></td><td>0</td><td>0</td><td>38 × 10<sup>-6</sup></td><td>6 × 10<sup>-6</sup></td><td>4 × 10<sup>-6</sup></td></tr><tr><td>Grain size (diff)</td><td>d</td><td>m</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td><td>0.001</td></tr><tr><td>Grain size exponent (diff)</td><td>m</td><td>–</td><td>2.0</td><td>2.0</td><td>3.0</td><td>0</td><td>0</td></tr></tbody></table>

horizontal advection, and vertical uplift (where in the latter two cases, the X and Y velocities and the Z velocities from ASPECT are used as model input) (58–60). Topography changes in the continental domain ($h \ge h_{sea}$) are therefore described by the following equation

$$ \frac{dh}{dt} = U - K_f A^m S^n + \frac{G}{A} \int_A \left( U - \frac{dh}{dt} \right) dA + K_c \nabla^2 h + \mathbf{v} \cdot \nabla h \quad (5) $$

where $h$ is the topographic elevation, $U$ the uplift rate, $K_f$ the bedrock erodibility, $A$ the drainage area, $S$ the slope, $m$ the drainage area exponent, $n$ the slope exponent, $G$ the deposition coefficient, $K_c$ the continental diffusion coefficient, and $\mathbf{v}$ the horizontal velocity. See Table 2 for the used parameter values.

Moreover, marine processes ($h < h_{sea}$) are described as follows

$$ \frac{dh}{dt} = K_m \nabla^2 h + Q_s + \mathbf{v} \cdot \nabla h + Q_o \quad (6) $$

<footnote>*Model input densities are scaled so that at surface temperatures ($T_0$ = 273 K or 0°C) these values are reached. diff = diffusion creep.</footnote>

<footnote>†Creep properties: dis = dislocation creep,</footnote>
<table><caption>Table 2. FastScape parameters.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Drainage area exponent</td><td>m</td><td>–</td><td>0.4</td></tr><tr><td>Slope exponent</td><td>n</td><td>–</td><td>1</td></tr><tr><td>Bedrock/sediment diffusivity</td><td>K<sub>c</sub></td><td>m<sup>2</sup> year<sup>-1</sup></td><td>5 &times; 10<sup>-3</sup></td></tr><tr><td>Bedrock/sediment erodibility</td><td>K<sub>f</sub></td><td>m<sup>0.2</sup> year<sup>-1</sup></td><td>1 &times; 10<sup>-5</sup></td></tr><tr><td>Bedrock/sediment deposition coefficient</td><td>G</td><td>–</td><td>1</td></tr><tr><td>Marine diffusivity</td><td>K<sub>m</sub></td><td>m<sup>2</sup> year<sup>-1</sup></td><td>200</td></tr><tr><td>Sand/shale ratio</td><td>F</td><td>–</td><td>1</td></tr><tr><td>Sand/shale porosity</td><td>&phi;</td><td>–</td><td>0</td></tr><tr><td>Sand/shale e-folding depth</td><td>z</td><td>m</td><td>0</td></tr><tr><td>Depth averaging thickness</td><td>L</td><td>m</td><td>100</td></tr><tr><td>Background sedimentation rate</td><td>Q<sub>0</sub></td><td>m year<sup>-1</sup></td><td>0</td></tr></tbody></table>

where $K_m$ represents the marine diffusion coefficient, $Q_s$ the sediment flux from the continent to the marine domain at the continent-marine boundary (60), and $Q_0$ a homogeneous marine background sedimentation rate.

## Geodynamic model design

### Model geometry

Our 2D rectangular domain is 450-km long and 200-km high and contains three layers together representing a 115-km-thick stable continental lithosphere: (i) a 20-km-thick wet quartzite upper crust (62), (ii) a 15-km-thick wet anorthite lower crust (63), and (iii) an 80-km-thick lithospheric mantle consisting of dry olivine (64, 65), overlying a wet olivine asthenosphere (66) (Fig. 7A and Table 1). Thicker upper crust (25 km) in the model domain center (between $x = 150$ km and $x = 300$ km, following a Gaussian distribution) produces a weaker strength profile and localizes kinematically driven tectonic deformation in this region (Fig. 7B). Furthermore, a randomized initial plastic (brittle) strain pattern with a maximum value of 0.5 is introduced in the upper 50 km of the crust to assist strain localization (Fig. 7A). The initial temperature follows a 1D steady-state continental geotherm in the lithosphere (67) and an adiabatic profile below.

### Viscoplastic rheology and weakening processes

We use a visco-plastic rheology (68) that combines diffusion creep, dislocation creep and Drucker-Prager plasticity. The friction angle is linearly weakened up to four times between brittle strain values of 0 and 1 (Fig. 7C and Table 1). Similarly, we apply a four times predicted viscous weakening between these same values of the viscous strain. Picard iterations on the Stokes solution are used to solve up to a nonlinear tolerance of $2 \times 10^{-5}$ for a maximum of 120 iterations per time step.

### Boundary conditions

A fixed temperature of 0°C is prescribed at the surface and a temperature of 1361°C at the bottom boundary. Tectonic deformation (divergence and convergence) is induced by prescribing outward and inward flow on the lateral model boundaries, respectively (Fig. 7A). This prescribed horizontal boundary velocity is $5 \text{ mm year}^{-1}$ on either lateral boundary, making for a total velocity of $10 \text{ mm year}^{-1}$

during either rifting or inversion. Flow through the lateral boundaries is compensated by inflow and outflow through the bottom boundary, which is controlled by a traction boundary condition set to a lithostatic pressure of $6.1 \times 10^9$ Pa (Fig. 7A). The latter boundary condition is imposed to conserve mass within the ASPECT model over time and to ensure isostatically balanced topography. Tectonic quiescence (leading to post-rift cooling) between divergence and convergence phases is simulated by prescribing a velocity of zero at the lateral model boundaries.

### FastScape parameters

We adopt a sea-level of 500 m below the initial surface level of the ASPECT model domain, a value that corresponds to the average elevation of stable continents with ~35-km crust far away from plume-related hotspot swells (69). All FastScape parameters are specified in Table 2. In particular, we apply a fluvial erosion coefficient ($K_f$) of $10^{-5} m^{0.2} year^{-1}$ and a continental hillslope diffusivity coefficient ($K_c$) of $5 \times 10^{-3} m^2 year^{-1}$ that are in agreement with previous model and observational constraints [see (59) and references therein].

### Discretization

The ASPECT mesh cell size is 10 km at the bottom of the domain, with a three-step refinement down to 1250 m above the 750°C isotherm, to properly capture brittle deformation in the colder parts of the lithosphere. We use second-order (Q2Q1) elements, and to optimally visualize the solution, results are output and analyzed on a grid with double the resolution of each cell (Fig. 7), the maximum resolution being 625 m.

### Varied model parameters

We systematically test the impact of plate tectonic parameters in our rift-inversion orogen models (Table 3). First, we apply different initial rifting durations of 5, 15, and 25 Myr, at a total plate motion velocity of $10 \text{ mm year}^{-1}$. These values roughly reflect the varying degrees of rifting before convergence [e.g., in the High Atlas, Pyrenees, and Alps, e.g., (22, 28, 70)], while the total convergence duration in all models is equal to the initial rifting duration plus 5 Myr to allow for initial collision to occur. Second, we apply a post-rift tectonic quiescence period of 0 or 20 Myr, which also reflect natural rift-inversion orogen settings (22, 28, 70).
<table><caption>Table 3. Variable model parameters.</caption><thead><tr><th>Model name (main text)</th><th>Model name [Materials and Methods section and supplementary data (37, 61)]</th><th>Rifting duration</th><th>Post-rift cooling duration</th><th>Inversion duration</th><th>Total model duration</th></tr></thead><tbody><tr><td>A</td><td>M1</td><td>5 Myr</td><td>0 Myr</td><td>10 Myr</td><td>15 Myr</td></tr><tr><td>-</td><td>M2</td><td>15 Myr</td><td>0 Myr</td><td>20 Myr</td><td>35 Myr</td></tr><tr><td>-</td><td>M3</td><td>25 Myr</td><td>0 Myr</td><td>30 Myr</td><td>55 Myr</td></tr><tr><td>-</td><td>M4</td><td>5 Myr</td><td>20 Myr</td><td>10 Myr</td><td>35 Myr</td></tr><tr><td>B</td><td>M5</td><td>15 Myr</td><td>20 Myr</td><td>20 Myr</td><td>55 Myr</td></tr><tr><td>C</td><td>M6</td><td>25 Myr</td><td>20 Myr</td><td>30 Myr</td><td>75 Myr</td></tr></tbody></table>

## Model analysis

Using the open source ParaView visualization software (www.paraview.org), we extract for each model the following three overall metrics, for each model data output interval of 500 kyr: the area of exhumed mantle (i.e., the total area of mantle above the initial Moho at 35 km depth; Figs. 1 and 7), the area of actively rising exhumed mantle (i.e., the area of exhumed mantle with upward vertical velocity > 0 m s⁻¹), and the area of mantle material within the adopted serpentinization window (i.e., the 200° to 350°C temperature window, in which serpentinization is considered to be most efficient if the water required for the reaction would be available).

Note that this serpentinization window extent is rather conservative as efficient serpentinization can also occur at higher temperatures and pressures, particularly during inversion and subduction [e.g., (19, 23)], but we base our analysis of natural H₂ potential on the temperature-dependent serpentinization formula from Liu et al. (12, 18). According to this formula, peak serpentinization efficiency (and thus bulk serpentinization) follows a bell curve that is centered in the 200° to 350°C serpentinization window, with very minor serpentinization occurring outside this window (Fig. 1C). This serpentinization window is based on conditions in rift systems (12). Although efficient serpentinization may also occur at higher pressures and temperatures [e.g., (19, 23)], we applied the same equation for assessing natural H₂ potential throughout the evolution of our modeled rift-inversion orogens to ensure an internally consistent analysis. We believe that this is permissible because tests with a larger serpentinization window (200° to 500°C) indicate that the main insights from our analysis remain robust (for more details, see also the “Nuance to our model analysis” section on nuances to our approach).

We subsequently adopt the general approach of Liu et al. (12) to calculate serpentinization capacity in kg(km year)⁻¹ and associated H₂ generation potential in mol(km year)⁻¹ in geodynamic models (Fig. 1C). This approach is based on the first-order assumption that active and sufficiently mature faults feature a high permeability, which allows for ample water circulation to enable high-degree serpentinization of mantle rocks at any moment in time. Therefore, at each output time step, we select the model mesh cells with mantle rocks that have accumulated a total plastic (i.e., brittle) strain ε ≥ 0.5 and that are undergoing considerable active deformation (second strain rate invariant ᾵ ≥ 2.5 × 10⁻¹⁵). Note that we also test different strain and strain rate thresholds (ε ≥ 0.1 or 1, and ᾵ ≥ 1.5 × 10⁻¹⁵ or 5 × 10⁻¹⁵) to estimate the variability in the results of our analysis

[see (37)]. We find that these reasonably different thresholds do not meaningfully change our conclusions.

Taking the selected model mesh cells with actively deforming mantle rocks, we apply the temperature-dependent serpentinization reaction formulas from (12, 18) to compute serpentinization rates in these mesh cells using ParaView's Python calculator application

$$ \frac{\partial \mathrm{Dserp}}{\partial t} = A e^{-[b_s(T - c_s)]^2} \qquad (7) $$

where $\frac{\partial \mathrm{Dserp}}{\partial t}$ is the serpentinization rate of a given model mesh cell (unit: s⁻¹) and T is the temperature in degrees Celsius of that cell, while empirically fitted kinetic coefficients ($A = 10^{-10}$ s⁻¹, $b_s = 2.5 \times 10^{-4}$°C⁻², $c_s = 270$°C) have been derived from experimental and theoretical observations on serpentinization kinematics (18).

With the serpentinization rate, we can calculate the serpentinization capacity ($M_t$) along the modelled section using the following equation from (12) in the ParaView Python calculator

$$ M_t = \sum_{i=1}^{i=n_{\mathrm{cell}}} \left( \frac{\partial \mathrm{Dserp}}{\partial t} \right) (i) S_i \rho_m \qquad (8) $$

where $(\frac{\partial \mathrm{Dserp}}{\partial t})(i)$ is the serpentinization rate of a given model mesh cell $i$, with $S_i$ being the area of that cell, and $\rho_m$ the density of the material in that cell. By summing over all cells (ncell), and given that our geodynamic ASPECT model is 2D, we obtain $M_t$ in units of kg(m s⁻¹), which we convert to kg (km year)⁻¹. Here, the “km” indicates length along-strike of the tectonic system (i.e., perpendicular to the 2D geodynamic ASPECT model; Figs. 2 and 7).

We subsequently derive the natural H₂ capacity in units of mol(km year)⁻¹, assuming that 1 kg of serpentinized mantle material generates between 100 and 300 mmol of natural H₂ (11). We thus obtain an upper limit estimate of the natural H₂ capacity in our modelled rift-inversion orogens while accounting for the considerable variations in natural H₂ generation potential that may occur in natural systems (11).

## Nuances to our model analysis

We must consider some nuances to our first-order modeling approach, which involves a number of assumptions and simplifications required for the models to be realized, as specified below. We believe that our results and the orogenic natural H₂ system concept (Fig. 5) we present in this study are robust, providing key insights and a
solid framework for understanding natural H₂ potential in rift-inversion orogens. Still, they should not be directly applied to interpret complex natural cases but rather serve as inspiration for more detailed assessments while keeping in mind the below nuances.

**Large-scale limitations to the geodynamic models**

Our models are essentially 2D (Fig. 2), whereas natural orogens are complex 3D structures [e.g., (71)], each with a unique history that requires assessment of its natural H₂ potential on a case-by-case basis, with an eye on variations along the strike of the system. Even so, the general structures we obtain in our models are realistic on a first-order basis (Figs. 1, 2, and 6). Moreover, because of constraints in the ASPECT code, melting and oceanic crust formation are not included in our models. However, we do develop the thermal profile that would be expected from an oceanic basin [Fig. 2, C and G, and (37)], whereas the thin sediment cover developing in rift basins provides a brittle layer that behaves in a similar way to oceanic crust. In addition, the mafic oceanic crust can also serpentinize, in a similar way as the mantle material exhumed in our modeled rift basins. Another caveat is that our analysis does not allow the detailed tracing of serpentinized and unserpentinized (fresh) mantle rock, so that some overlap and thus overestimation of serpentinized mantle volumes may occur. Still, our models suggest general exhumation of fresh mantle material during rift-inversion orogen evolution: by mantle uprising below the rift basin during rifting and by overthrusting of the overriding mantle wedge during inversion (Fig. 3).

**Extent of the serpentinization window and omittance of subduction-related serpentinization**

For our natural H₂ potential analysis, we focus on serpentinization of exhumed mantle material, which can be considered to be most efficient at temperatures between 200° and 350°C during rifting (12, 18). To ensure consistency in our analysis, we apply this same window to assess the natural H₂ potential throughout the whole evolution of the rift-inversion orogens. Yet, this serpentinization window is rather conservative when considering the inversion stage of rift-inversion evolution because efficient serpentinization can also occur at higher temperatures and pressures (19, 23). We conducted a test with a more extensive serpentinization window of 200° to 500°C that may be more representative for deeper (10 km) serpentinization during inversion (19). We find that, although there would be more serpentinization in the system overall, the ratio of serpentinization per time unit between rifting and inversion would remain similar or would even be larger when applying these different windows for rifting and inversion, respectively [see details in (37)]. Regardless of how the serpentinization window is parameterized, the fact that orogens are colder than rifts means that the orogenic serpentinization window is decisively larger and hence that considerably more natural H₂ is generated.

Moreover, because we focus our analysis on exhumed mantle material (i.e., in the overriding mantle wedge in a rift-inversion orogen; Fig. 1C), we omit the potential additional serpentinization and natural H₂ generation in the deepest parts of the system during inversion (at depths of 50 to 60 km in the subduction channel) (19, 23). Including this deep inversion-related serpentinization, a challenge beyond the scope of our study, would increase the already distinct difference in serpentinization between rifting and inversion.

**Water availability**

Liquid water is key to natural H₂ generation because serpentinization cannot occur without it. In our analysis, we assume that liquid

water, either from surface (marine or meteoric) or from deep sources (metamorphic), is readily available along large active faults that cross exhumed mantle material. This assumption is reasonable during rifting, when seawater can percolate along large normal faults. The situation is more complex in rift-inversion orogens, where meteoric water may reach deep rocks as a result of the high hydraulic head and temperature differences (21), whereas dehydration reactions linked to subduction may release large volumes of water that can induce serpentinization in the overlying mantle rocks as well [e.g., (23) (Fig. 5)]. A thorough understanding of fluid flow and water-rock reactions is critical for understanding the spatial and temporal distribution of serpentinization in natural systems. For example, the impact of density changes during serpentinization requires attention because the volumetric expansion of serpentinizing mantle material may either clog water pathways (72) or induce new fracturing and increase permeability (73). Still, we can expect large-scale faults to generally provide efficient pathways for water circulation, both in rifting and inversion settings.

**Quantities of generated natural H₂**

Serpentinization cannot directly be equated to H₂ generation since next to temperature, pressure and water availability, the type (composition) and freshness of mantle material, as well as the origin and associated chemistry of the water involved in serpentinization modulate natural H₂ generation potential (11, 40, 68). As a result, serpentinization reactions are very complex, and corresponding natural H₂ generation from mantle material has a considerable range in nature. We account for this uncertainty by applying a range of 100 to 300 mmol of natural H₂ per kilogram of mantle rock in our first-order natural H₂ potential analysis (11).

**Natural H₂ migration and preservation**

Other limitations to our models concern the omission of detailed sedimentary history (i.e., timing and availability of reservoirs and seals) and of detailed natural H₂ migration pathways to potential reservoirs because highest model resolution of 625 m only allows the tracing of the largest faults [e.g., (21, 25)]. Rock-water interactions and the aforementioned volumetric expansion of serpentinizing mantle material and impact on fault permeability (either increasing or decreasing) also affect the efficiency of natural H₂ migration pathways (72, 73), although the size of H₂ molecules means that they can move more easily along tight passages than other substances. These reactions are also relevant for the quality of reservoir rocks and seals. A further important question is whether the high buoyancy of H₂ molecules allows them to move upward with relative freedom, or whether their migration will be meaningfully affected and redirected by fluid flows along major faults.

Moreover, H₂ gas is highly reactive and may easily be consumed by (bio)chemical reactions during serpentinization, as well as on its way to, or within reservoirs [e.g., (8, 74)]. For example, during or shortly after serpentinization, the resulting natural H₂ can be swiftly converted to CH₄ or other by-products (23, 51), which may subsequently end up in the reservoirs that would otherwise accumulate the natural H₂. Last, although we point out that reservoir temperatures of natural H₂ fields should ideally be higher than the 122°C biotic fringe at several kilometer depth for the best preservation conditions (20, 40), much shallower natural H₂ accumulations at some hundreds of meters depth are known from Mali and Albania (16, 33), potentially due to very high H₂ fluxes that replenish any loss of H₂.
REFERENCES AND NOTES

1. L. Truche, E. F. Bazarkina, Natural hydrogen the fuel of the 21st century. *E3S Web Conf. 98* 03006 (2019).

2. R. W. Howarth, M. Z. Jacobson, How green is blue hydrogen? *Energy Sci. Eng.* 9, 1676–1687 (2021).

3. B. S. Lollar, T. C. Onstott, G. Lacrampe-Couloume, C. J. Ballentine, The contribution of the Precambrian continental lithosphere to global H₂ production. *Nature* **516**, 379–382 (2014).

4. V. Zgonnik, The occurrence and geoscience of natural hydrogen: A comprehensive review. *Earth Sci. Rev.* **203**, 103140 (2020).

5. N. J. P. Smith, T. J. Shepherd, M. T. Styles, G. M. Williams, Hydrogen exploration: A review of global hydrogen accumulations and implications for prospective areas in NW Europe. *Geol. Soc. London Petrol. Geol. Conf. Ser.* **6**, 349–358 (2005).

6. E. C. Gaucher, New perspectives in the industrial exploration for native hydrogen. *Elements* **16**, 8–9 (2020).

7. E. C. Gaucher, I. Moretti, N. Pelissier, G. Burridge, N. Gonthier, The place of natural hydrogen in the energy transition: A position paper. *Eur. Geol. J.* **55**, 5–9 (2023).

8. F. Klein, J. D. Tarnas, W. Bach, Abiotic sources of molecular hydrogen on Earth. *Elements* **16**, 19–24 (2020).

9. R. W. Wilson, G. A. Houseman, S. J. H. Buiter, K. J. W. McCaffrey, A. G. Doré, Fifty years of the Wilson cycle concept in plate tectonics: An overview. *Geol. Soc. Lond. Spec. Publ.* **470**, 1–17 (2019).

10. A. S. Merdith, P. del Real, I. Daniel, M. Andreani, N. M. Wright, N. Cottice, Pulsated global hydrogen and methane flux at mid-ocean ridges driven by Pangea breakup. *Geochem. Geophys. Geosys.* **21**, e2019GC008869 (2020).

11. E. Albers, W. Bach, W. M. Pérez-Gussinyé, C. Cammon, T. Frederichs, Serpentinitization-driven H₂ production from continental break-up to mid-ocean ridge spreading: Unexpected high rates at the West Iberia margin. *Front. Earth Sci.* 9, 673063 (2021).

12. Z. Liu, M. Perez-Gussinyé, J. García-Pintado, L. Mezri, W. Bach, Mantle serpentinization and associated hydrogen flux at North Atlantic magma-poor rifted margins. *Geology* **51**, 284–289 (2023).

13. A. S. Merdith, I. Daniel, D. Sverjensky, M. Andreani, B. Mather, S. Williams, A. Vitale Brovarone, Global hydrogen production during high-pressure serpentinization of subducting slabs. *Geochem. Geophys. Geosyst.* **24**, e2023GC010947 (2023).

14. J. A. Leong, M. Nielsen, N. McQueen, R. Karolyté, D. J. Hillegonds, C. Ballentine, T. Darrah, W. McGillis, P. Kelemen, H₂ and CH₄ outgassing rates in the Samail ophiolite, Oman: Implications for low-temperature, continental serpentinization rates. *Geochem. Cosmochim. Acta* **347**, 1–15 (2023).

15. E. Deville, A. Prinzhofer, The origin of N₂-H₂-CH₄-rich natural gas seepages in ophiolitic context: A major and noble gases study of fluid seepages in New Caledonia. *Chem. Geol.* **440**, 139–147 (2016).

16. L. Truche, F.-V. Donzé, E. Goskolli, B. Muceku, C. Loisy, C. Monnin, H. Dutoit, A. Cerepi, A deep reservoir for hydrogen drives intense degassing in the Bulqizë ophiolite. *Science* **383**, 618–621 (2024).

17. D. Lévy, M. Boka-Mene, A. Meshi, I. Fejza, T. Guermont, B. Hauville, N. Pelissier, Looking for natural hydrogen in Albania and Kosova. *Front. Earth Sci.* **11**, 1167634 (2023).

18. S. Emmanuel, B. Berkowitz, Suppression and stimulation of seafloor hydrothermal convection by exothermic mineral hydration. *Earth Planet. Sci. Lett.* **243**, 657–668 (2006).

19. R. Huang, W. Sun, X. Ding, Z. Yusheng, M. Song, Effect of pressure on the kinetics of peridotite serpentinization. *Phys. Chem. Miner.* **47**, 33 (2020).

20. A. Vinsot, C. A. J. Appelo, M. Lundy, S. Wechner, Y. Lettry, C. Lerouge, A. M. Fernández, M. Labat, C. Tournassat, P. De Canniere, B. Schwyn, J. McKelvie, S. Dewonck, P. Bossart, J. Delay, In situ diffusion test of hydrogen gas in the Opalinus Clay. *Geol. Soc. Spec. Pub.* **400**, 563–578 (2014).

21. N. Lefeuvre, L. Truche, F.-V. Donzé, F. Gal, J. Tremosa, R.-A. Fakoury, S. Calassou, E. C. Gaucher, Natural hydrogen migration along thrust faults in foothill basins: The North Pyrenean Frontal Thrust case study. *Appl. Geochem. Geophys.* **145**, 105396 (2022).

22. S. Ferrando, M. L. Frezzotti, P. Orione, R. C. Conte, R. Compagnoni, Late-Alpine rodingitization in the Bellecombe meta-ophiolites (Aosta Valley, Italian Western Alps): Evidence from mineral assemblages and serpentinization-derived H₂-bearing brine. *Int. Geol. Rev.* **52**, 1220–1243 (2010).

23. A. Vitale Brovarone, D. A. Sverjensky, F. Piccoli, F. Ressico, D. Giovanelli, I. Daniel, Subduction hides high-pressure sources of energy that may feed the deep subsurface biosphere. *Nat. Commun.* **11**, 3880 (2020).

24. F. Giuntoli, L. Menegon, G. Siron, F. Cogni, H. Leroux, R. Compagnoni, M. Rossi, A. Vitale Brovarone, Methane-hydrogen-rich fluid migration may trigger seismic failure in subduction zones at forearc depths. *Nat. Commun.* **15**, 480 (2024).

25. N. Lefeuvre, L. Truche, F.-V. Donzé, M. Ducoux, G. Barré, R.-A. Fakoury, S. Calassou, E. C. Gaucher, Native H₂ exploration in the Western Pyrenean Foothills. *Geochem. Geophys. Geosys.* **22**, e2021GC009917 (2021).

26. F. Zwaan, P. Pilz, S. Niedermann, M. Zimmer, N. Lefeuvre, V. Petit, A. Vieth-Hillebrand, E. C. Gaucher, C. Schmidt-Hattenberger, S. Brune, Soil and spring water gas data from an area

with natural hydrogen (H₂) degassing, NW Pyreneaneland, France (GFZ Data Services, 2024).

27. D. A. Vasey, J. B. Naliboff, E. Cowgill, S. Brune, A. Glerum, F. Zwaan, Impact of rift history on the structural style of intracontinental rift-inversion orogens. *Geology* **52**, 429–434 (2024).

28. S. M. Schmid, E. Kissling, T. Diehl, D. J. J. van Hinsbergen, G. Molli, Ivrea mantle wedge, arc of the Western Alps, and kinematic evolution of the Alps–Apennines orogenic system. *Swiss J. Geosci.* **110**, 581–612 (2017).

29. M. Ford, E. Masini, J. Vergés, R. Pik, S. Ternois, J. Léger, A. Dielforder, G. Frasca, A. Grool, C. Vinciguerra, T. Bernard, P. Angrand, A. Crémades, G. Manatschal, S. Chevrot, L. Jolivet, F. Mouthereau, I. Thinon, S. Calassou, Evolution of a low convergence collisional orogen: A review of Pyrenean orogenesis. *Bull. Soc. Géol. Fr.* **193**, 19 (2022).

30. S. Jammes, R. S. Huismans, Structural styles of mountain building: Controls on lithospheric rheologic stratification and extensional inheritance. *J. Geophys. Res. Solid Earth* **117**, 10.1029/2012JB009376 (2012).

31. Z. Erdós, R. S. Huismans, P. van der Beek, C. Thieulot, Extensional inheritance and surface processes as controlling factors of mountain belt structure. *J. Geophys. Res. Solid Earth* **119**, 9042–9061 (2014).

32. A. Jourdon, L. Le Pourhiet, F. Mouthereau, E. Masini, Role of rift maturity on the architecture and shortening distribution in mountain belts. *Earth Planet. Sci. Lett.* **512**, 89–99 (2019).

33. A. Prinzhofer, C. S. T. Cissé, A. B. Diallo, Discovery of a large accumulation of natural hydrogen in Bourakebougou (Mali). *Int. J. Hydrogen Energy* **43**, 19315–19326 (2018).

34. M. Kronbichler, T. Heister, W. Bangerth, High accuracy mantle convection simulation through modern numerical methods. *Geophys. J. Int.* **191**, 12–29 (2012).

35. T. Heister, J. Dannberg, R. Gassmöller, W. Bangerth, High accuracy mantle convection simulation through modern numerical methods—II: Realistic models and problems. *Geophys. J. Int.* **210**, 833–851 (2017).

36. D. Neuharth, S. Brune, T. Wrona, A. Glerum, J. Braun, X. Yuan, X., Evolution of rift systems and their fault networks in response to surface processes. *Tectonics* **41**, e2021TC00716 (2022).

37. F. Zwaan, S. Brune, A. C. Glerum, D. A. Vasey, J. B. Naliboff, G. Manatschal, E. C. Gaucher, Results of serpentinization and associated natural H₂ potential analysis in rift-inversion orogen systems, using a numerical geodynamic modelling approach (GFZ Data Services, 2025). https://doi.org/10.5880/GFZ.2.5.2025.001.

38. M. de Saint Blanquat, F. Bajolet, A. Grand’Homme, A. Proietti, M. Zanti, A. Boutin, C. Clerc, Y. Lagabrielle, P. Labaume, Cretaceous mantle exhumation in the central Pyrenees: New constraints from the peridotites in eastern Ariège (North Pyrenean zone, France). *C. R. Geosci.* **348**, 268–278 (2016).

39. C. Tichadou, M. Godard, M. Muñoz, P. Labaume, A. Vauchez, E. C. Gaucher, S. Calassou, Mineralogical and geochemical study of serpentinized peridotites from the North-Western Pyrenees: New insights on serpentinization along magma-poor continental passive margins. *Lithos* **406–407**, 106521 (2021).

40. K. Takei, K. Nakamura, T. Toki, U. Tsunogai, M. Miyazaki, J. Miyazaki, H. Hirayama, S. Nakagawa, T. Nunoura, K. Horikoshi, Cell proliferation at 122°C and isotopically heavy CH₄ production by a hyperthermophilic methanogen under high-pressure cultivation. *Proc. Natl. Acad. Sci. U.S.A.* **105**, 10949–10954 (2008).

41. F. Osselin, C. Soulaine, C. Fauguerolles, E. C. Gaucher, B. Scaillet, M. Pichavant, Orange hydrogen is the new green. *Nat. Geosci.* **15**, 765–769 (2022).

42. M. W. Downey, Oil and natural gas exploration, in *Encyclopedia of Energy* (Elsevier, 2004), pp. 549–558.

43. J. A. Muñoz, Alpine Orogeny: Deformation and structure in the Northern Iberian Margin (Pyrenees s.l.), in *The Geology of Iberia: A Geodynamic Approach: Volume 3: The Alpine Cycle, Regional Geology Studies* (Springer International Publishing, 2019), pp. 433–451.

44. Y. Wang, S. Chevrot, V. Monteiller, D. Komatitsch, F. Mouthereau, G. Manatschal, M. Sylvander, J. Diaz, M. Ruiz, F. Grimaud, S. Benahmed, H. Pauchet, R. Martin, The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic P waves. *Geology* **44**, 475–478 (2016).

45. M. Lehujeur, S. Chevrot, A. Villaseñor, E. Masini, N. Saspiturry, E. Lescoutre, M. Sylvander, Three-dimensional shear velocity structure of the Mauléon and Arzacq Basins (Western Pyrenees). *Bull. Soc. Géol. Fr.* **192**, 47 (2021).

46. S. Chevrot, M. Sylvander, A. Villaseñor, J. Diaz, L. Stehly, P. Boulé, V. Monteiller, R. Martin, M. Lehujeur, S. Beller, J. Bribes, A. Bitri, S. Calassou, M. Collin, M. Ford, L. Jolivet, G. Manatschal, E. Masini, F. Mouthereau, O. Vidal, Passive imaging of collisional orogens: a review of a decade of geophysical studies in the Pyrénées. *Bull. Soc. Géol. Fr.* **193**, 1 (2022).

47. "Global Change Data Lab," Our World in Data (2024); https://ourworldindata.org/energy#explore-data-on-energy.

48. S. M. Schmid, D. Bernoulli, F. Bügenschuh, L. Matenco, S. Schefer, R. Schuster, M. Tischler, K. Ustaszewski, The Alpine-Carpathian-Dinaridic orogenic system: Correlation and evolution of tectonic units. *Swiss J. Geosci.* **101**, 139–183 (2008).

49. A. H. F. Robertson, B. Trivić, N. Đerić, I. Bucur, Tectonic development of the Vardar ocean and its margins: Evidence from the republic of Macedonia and Greek Macedonia. *Tectonophysics* **595–596**, 25–54 (2013).
50. S. Mazzoli, A. M. Algarra, Deformation partitioning during transpressional emplacement of a 'mantle extrusion wedge': The Ronda peridotites, western Betic Cordillera, Spain. *J. Geol. Soc.* **168**, 373–382 (2011).

51. G. Etiope, I. Vadillo, M. J. Whiticar, M. Marques, P. M. Carreira, I. Tiago, J. Benavente, P. Jiménez, B. Urresti, Abiotic methane seepage in the Ronda peridotite massif, southern Spain. *Appl. Geochem.* **66**, 101–113 (2016).

52. H. Furnes, Y. Dilek, G. Zhao, I. Safonova, M. Santosh, Geochemical characterization of ophiolites in the Alpine-Himalayan Orogenic Belt: Magmatically and tectonically diverse evolution of the Mesozoic Neotethyan oceanic crust. *Earth Sci. Rev.* **208**, 103258 (2020).

53. I. Rose, B. Buffett, T. Heister, Stability and accuracy of free surface time integration in viscous flows. *Phys. Earth Planet. Int.* **262**, 90–100 (2017).

54. R. Gassmöller, H. Lokavarapu, E. Heien, E. G. Puckett, W. Bangerth, Flexible and scalable particle-in-cell methods with adaptive mesh refinement for geodynamic computations. *Geochem. Geophys. Geosys.* **19**, 3596–3604 (2018).

55. A. Glerum, C. Thieulot, M. Fraters, C. Blom, W. Spakman, Nonlinear viscoplasticity in viscoplasticity in ASPECT: Benchmarking and applications to subduction. *Solid Earth* **9**, 267–294 (2018).

56. W. Bangerth, J. Dannberg, M. Fraters, R. Gassmoeller, A. Glerum, T. Heister, R. Myhill, J. Naliboff, ASPECT v2.4.0 (v2.4.0) (Zenodo, 2022).

57. W. Bangerth, J. Dannberg, M. Fraters, R. Gassmoeller, A. Glerum, T. Heister, R. Myhill, J. Naliboff, ASPECT: Advanced solver for problems in Earth's Convection user manual. (Figshare; 2023); https://doi.org/10.6084/M9.FIGSHARE.4865333.

58. J. Braun, S. D. Willett, A very efficient $O(n)$, implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution. *Geomorphology* **180–181**, 170–179 (2013).

59. X. P. Yuan, J. Braun, L. Guerit, D. Rouby, G. Cordonnier, A new efficient method to solve the stream power law model taking into account sediment deposition. *J. Geophys. Res. Earth* **124**, 1346–1365 (2019).

60. X. P. Yuan, J. Braun, L. Guerit, B. Simon, B. Bovy, D. Rouby, D. C. Robin, R. Jiao, Linking continental erosion to marine sediment transport and deposition: A new implicit and $O(N)$ method for inverse analysis. *Earth Planet. Sci. Lett.* **524**, 115728 (2019).

61. F. Zwaan, S. Brune, A. C. Glerum, D. A. Vasey, J. B. Naliboff, G. Manatschal, E. C. Gaucher, Rift-inversion orogen modelling files for natural $H_2$ potential assessment (Zenodo, 2025); https://doi.org/10.5281/zenodo.14672886.

62. E. H. Rutter, K. H. Brodie, Experimental grain size-sensitive flow of hot-pressed Brazilian quartz aggregates. *J. Struct. Geol.* **26**, 2011–2023 (2004).

63. E. Rybacki, M. Gottschalk, R. Wirth, G. Dresen, Influence of water fugacity and activation volume on the flow properties of fine-grained anorthite aggregates. *J. Geophys. Res.* **111**, B03203 (2006).

64. I. M. Artemieva, Global 1° × 1° thermal model TC1 for the continental lithosphere: Implications for lithosphere secular evolution. *Tectonophysics* **416**, 245–277 (2006).

65. M. E. Pasyanos, T. G. Masters, G. Laske, Z. Ma, Z., LITH1.0: An updated crust and lithospheric model of the Earth. *J. Geol. Res. Solid. Earth* **119**, 2153–2173 (2014).

66. G. Hirth, D. Kohlstedt, Rheology of the upper mantle and the mantle wedge: A view from the experimentalists, in *Inside the subduction factory*, J. Eiler, Ed. (Blackwell Publishing Ltd, 2004).

67. D. S. Chapman, Thermal gradients in the continental crust. *Geol. Soc. London Spec. Publ.* **24**, 63–70 (1986).

68. H. M. Lamadrid, J. D. Rimstidt, E. M. Schwarzenbach, F. Klein, S. Ulrich, A. Dolocan, R. J. Bodnar, Effect of water activity on rates of serpentinization of olivine. *Nat. Commun.* **8**, 16107 (2017).

69. T. Theunissen, R. S. Huismans, G. Lu, N. Riel, Relative continent–mid-ocean ridge elevation: A reference case for isostasy in geodynamics. *Earth Sci. Rev.* **233**, 104153 (2022).

70. W. Beauchamps, R. W. Allmendinger, M. Barazangi, A. Demnati, M. El Alji, M. Dahmani, Inversion tectonics and the evolution of the high Atlas Mountains, Morocco, based on a geological-geophysical transect. *Tectonics* **18**, 163–184 (1999).

71. G. Manatschal, P. Chenin, R. Lescoutre, J. Miró, P. Cadenas, N. Saspiturry, E. Masini, S. Chevrot, M. Ford, L. Jolivet, F. Mouthereau, I. Thinon, B. Issautier, S. Calassou, The role of

inheritance in forming rifts and rifted margins and building collisional orogens: A Biscay-Pyrenean perspective. *Bull. Soc. Géol. Fr.* **192**, 55 (2021).

72. A. Farough, D. E. Moore, D. A. Lockner, R. P. Lowell, Evolution of fracture permeability of ultramatic rocks undergoing serpentinization at hydrothermal conditions: An experimental study. *Geochem. Geophys. Geosys.* **17**, 44–55 (2016).

73. M. Uno, K. Koyanagawa, H. Kasahara, A. Okamoto, N. Tsuchiya, Volatile-consuming reactions fracture rocks and self-accelerate fluid flow in the lithosphere. *Proc. Natl. Acad. Sci. U.S.A.* **119**, e2110776118 (2022).

74. S. L. Worman, L. F. Pratson, J. A. Karson, W. H. Schlesinger, Abiotic hydrogen ($H_2$) sources and sinks near the mid-ocean ridge (MOR) with implications for the subseafloor bio-sphere. *Proc. Natl. Acad. Sci. U.S.A.* **117**, 13283–13293 (2020).

75. F. Kolawole, J. C. Evenick, Global distribution of geothermal gradients in sedimentary basins. *Geosci. Front.* **14**, 101685 (2023).

76. Woudloper, Tectonic map Mediterranean EN (Wikimedia Commons, 2005); https://en.m.wikipedia.org/wiki/File:Tectonic_map_Mediterranean_EN.shtml

77. C. Faccenna, T. W. Becker, L. Auer, A. Billi, L. Boschi, J. P. Brun, F. A. Capitanio, F. Funiciello, F. Horváth, L. Jolivet, C. Piromallo, L. Royden, F. Rossetti, E. Serpelloni, Mantle dynamics in the Mediterranean. *Rev. Geophys.* **52**, 283–332 (2014).

78. J.A. Muñoz, The Pyrenees, in *The Geology of Spain*, W. Gibbons, T. Moreno, Eds. (Geological Society of London, 2002), pp. 370–385.

79. M. L. Arboleya, A. Teixell, M. Charroud, M. Julivert, A structural transect through the high and middle Atlas of Morocco. *J. Afr. Earth Sci.* **39**, 319–327 (2004).

80. C. L. Rosenberg, E. Kissling, Three-dimensional insight into Central-Alpine collision: Lower-plate or upper-plate indentation? *Geology* **41**, 1219–1222 (2013).

**Acknowledgments:** We thank our colleagues at the GFZ Helmholtz Centre for Geosciences in Potsdam who provided valuable support and discussion. We also thank the Computational Infrastructure for Geodynamics (https://geodynamics.org), funded by the National Science Foundation under award EAR-0949446 and EAR-1550901, for supporting the development of ASPECT. We gratefully acknowledge L. Truche, the reviewers for their helpful feedback on our manuscript, as well as the editor for organizing the review process. **Funding:** This research was principally funded through a GFZ Discovery Fellowship awarded to F.Z. S.B. was funded by the European Union (ERC, EMERGE, 101087245). A.G. received support from the Helmholtz Association through the Helmholtz Recruitment Initiative 0316. G.M. was supported by the Interdisciplinary Thematic Institute Geosciences for the energy system transition, as part of the ITI 2021-2028 program of the University of Strasbourg, CNRS and Inserm, funded by IdEx Unistra (ANR-10-IDEX-0002), and by SFRI-STRAT'US project (ANR-20-SFRI-0012). E.G.'s working time was covered by Lavoisier $H_2$ Geoconsult's internal HNAt project. Computing time was granted by the Resource Allocation Board and provided on the supercomputer Lise and Emmy at NHR@ZIB and NHR@Göttingen as part of the NHR infrastructure. The calculations for this research were conducted with computing resources under project no. bbp00039. **Author contributions:** Conceptualization: F.Z., S.B., G.M., and E.G. Methodology: F.Z., A.G., S.B., D.V., and J.N. Software: F.Z., A.G., S.B., D.V., and J.N. Validation: F.Z. and S.B. Formal analysis: F.Z. and S.B. Investigation: F.Z. Resources: S.B. and A.G. Data Curation: F.Z. and A.G. Writing—original draft: All authors. Writing—review and editing: All authors. Visualization: F.Z. Project administration: F.Z. and S.B. Funding acquisition: F.Z., S.B., G.M., and E.G. **Competing interests:** The authors declare that they have no competing interests. **Data and materials availability:** Supplementary material containing modeling software, input files, and model analysis scripts is publicly available at https://doi.org/10.5281/zenodo.14672886 (61). Moreover, additional Supplementary Material, in the form of a GFZ data publication containing additional details on model results as well as videos of each model, is publicly available at https://doi.org/10.5880/GFZ.2.5.2025.001 (37). All other data needed to evaluate the conclusions of our study are present in the paper.

Submitted 2 July 2024

Accepted 17 January 2025

Published 19 February 2025

10.1126/sciadv.adr3418