# Water budget in flat vs. steep subduction: implication for volcanism and potential for H2 production

Alexis Gauthier, Tiphaine Larvet, Laetitia Le Pourhiet, Isabelle Moretti

► To cite this version:

Alexis Gauthier, Tiphaine Larvet, Laetitia Le Pourhiet, Isabelle Moretti. Water budget in flat vs. steep subduction: implication for volcanism and potential for H2 production. Bulletin de la Société Géologique de France, 2024, 195, pp.26. 10.1051/bsgf/2024026. hal-04880223

HAL Id: hal-04880223

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

Submitted on 10 Jan 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>

Distributed under a Creative Commons Attribution 4.0 International License
# Water budget in flat vs. steep subduction: implication for volcanism and potential for H₂ production

Alexis Gauthier<fnref n="1" />,*, Tiphaine Larvet<fnref n="1" /> ID, Laetitia Le Pourhiet<fnref n="1" />,<fnref n="2" /> ID and Isabelle Moretti<fnref n="1" />,<fnref n="3" /> ID

Received: 29 January 2024 / Accepted: 21 October 2024 / Publishing online: 23 December 2024

**Abstract** – In this study, we model the processes of (de)hydration and melting within subduction zones using a thermo-mechanical modeling approach. Multiple 2D simulations are conducted to investigate how the subduction angle influences the water budget in oceanic-continental subduction, focusing mainly on the variation of slab dip angle along the strike of the Andes. It appears that in the case of flat subduction, the mantle hydration zone is large, extending up to 500 km from the trench. This extent depends on the length of flat slab segment which, in turn, depends on the velocity of the overriding plate. In the case of a steep subduction, the zone is narrower, and is located between the trench and the volcanic arc. Magma formation competes with hydration of the mantle wedge for the use of water expelled from the subducting plate. In the transition from a steep to a flat slab, the mantle hydration zone widens and the volcanic zone moves away from the trench. The oceanic crust may undergo melting, leading to a change in magma composition and the development of adakitic volcanism, before volcanism diminishes in intensity and then disappears. Our study provides geodynamic insights into observations related to volcanism in the Pampean flat slab in South America. Using the quantification of water involved in mantle wedge hydration as a proxy for H₂ production, we propose that flat subductions are the most promising areas for H₂ exploration. Additionally, deep H₂ production appears to be particularly sensitive to the amount of subducted sediments, regardless of whether subduction is flat or steep. Lower plate serpentinization does not affect deep H₂ production.

**Keywords:** Subduction / water budget / numerical modeling / serpentinization / Andes / magmatism / natural hydrogen

**Résumé** – Comparaison du bilan d'eau entre les subductions classiques et plates : implication pour le volcanisme et le potentiel pour la production d'H₂. Dans cette étude, nous modélisons les processus de (dé) hydration et de fusion au sein des zones de subduction en utilisant une approche de modélisation thermo- mécanique. Plusieurs simulations 2D sont menées pour étudier comment le pendage du panneau plongeant influence le bilan hydrique lors de la subduction océan-continent, dans le but de reproduire les variations du pendage le long des Andes. Il apparaît qu'en cas de subduction plate la zone d'hydratation du manteau est large et peut s'étendre jusqu'à 500 km de la fosse. L'étendue dépend de la longueur du segment de slab plat, qui est-elle- même influencée par la vitesse de la plaque chevauchante. En cas de slab penté, la zone est plus étroite, et elle est située entre la fosse et l'arc volcanique. La formation de magma entre en compétition avec l'hydratation du coin du manteau pour l'utilisation de l'eau expulsée de la plaque plongeante. Lors du passage d'un slab penté à un slab plat, la zone d'hydratation du manteau s'élargit et la zone volcanique s'éloigne de la fosse. La croûte océanique de la plaque plongeante peut fondre, entrainant une variation de la composition du magma avant que le volcanisme ne diminue en intensité puis ne disparaisse. Notre étude offre de nouvelles perspectives géodynamiques sur les observations relatives au volcanisme dans la région du Pampean flat slab en Amérique du Sud. En utilisant la quantification de l'eau impliquée dans l'hydratation du coin mantellique comme proxy potentiel pour la production d'H₂, nous proposons que les zones de subduction plate soient les plus prometteuses pour l'exploration d'H₂. De plus, il semble que la production profonde d'H₂ soit particulièrement sensible à la quantité de sédiments subductés, que la subduction soit plate ou abrupte. La serpentinisation de la plaque inférieure n'affecte pas la production profonde d'H₂.

**Mots clés** : subduction / bilan d'eau / modélisation numérique / serpentinisation / Andes / magmatisme / hydrogène naturel

<footnote n="1">¹ Sorbonne Université, Institut des Sciences de la Terre, ISTE, 75005 Paris, France</footnote>

<footnote n="2">² Institut Universitaire de France</footnote>

<footnote n="3">³ Université Pau Pays de l'Adour, UPPA, 64000 Pau, France</footnote>

<footnote>*Corresponding author: alexis.gauthier@sorbonne-universite.fr</footnote>
# 1 Introduction

## 1.1 Geodynamic context

The subduction between the Nazca Plate and South America is characterized by rapid convergence, driven by the expansion of both the Pacific and Atlantic oceans. This results in absolute motion of both oceanic and continental plate toward the trench (Sdrolias and Müller, 2006). Our study focuses on the segment between the paleo-flat slab of Puna and the current flat slab of the Pampean Ranges (Ramos and Folguera, 2009; Martinod et al., 2020) (20°S–30°S) over the past 20 million years.

The Pampean flat slab (Fig. 1), which began during early Miocene, around 18 Ma ago, with the reduction of the subduction angle of the Nazca Plate beneath South America (Litvak et al., 2007). Geological observations show that this reduction in slab dip angle does not mark the end of magmatic activity in the area but is associated with a transient period including (1) landward migration of the volcanic arc (Kay and Abbruzzi, 1996; Litvak et al., 2007; Ramos et al., 2002), (2) a decrease in the amount of lava expelled by volcanoes (Ramos and Folguera, 2009), (3) change in magma chemistry (increase in the La/Yb) with the presence of adakitic-type volcanoes (Kay and Mpodozis, 2002; Poma et al., 2017, 2023) until cessation 4.5 Ma ago (Poma et al., 2017). Adakitic volcanism can be initiated by various processes (Kay and Kay, 2002). Litvak et al. (2007) and Poma et al. (2017, 2023) suggested that the adakitic volcanism in the Pampean flat slab region is caused by magma contamination from the thickened continental crust. In contrast, Gutscher et al. (2000), using a thermodynamic approach, argue that the adakitic signature observed in volcanoes during the Pampean flat slab phase is the result of oceanic crust melting in a flat subduction environment.

Changes in volcanic activity are associated with flat slab segments, as previous studies have shown a correlation between these segments and the absence of active volcanoes across South America (Ramos and Folguera, 2009; Fig. 1). Similarly, Moretti et al. (2023) have detected a variation in the H₂ emanation function of the slab dip angle in the south Bolivian altiplano. While in Colombia (Fig. 1), accretion has been active and the presence of Mesozoic ultramafic rocks leads to a diversity of possible H₂ origin (Carrillo Ramirez et al., 2023), southward from Ecuador (Fig. 1), where there is no obducted mafic material, the variation of H₂ measured in the south Bolivian altiplano by Moretti et al. (2023) is likely to be derived from hydration of the mantle wedge. Since volcanic activity, and H₂ emanations alike, are closely linked to water circulations, we aim here at investigating how the water budget of flat subduction zones differs significantly from that of steep subduction zones.

## 1.2 Water in subduction zones

During oceanic-continental subduction, the downgoing plate releases water as it sinks into the Earth's mantle. For the oceanic plate, bounded water within the minerals is initially contained in the sediments, the crust and the first few kilometers of the lithospheric mantle (Faccenda, 2014). In sediments, water is mainly transported by lawsonite and

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Present flat slab segment and volcanism in South America, modified from (Ramos and Folguera, 2009).</figcaption></figure>

phengite at depth. Phengite remains stable up to 10 GPa (i.e. about 330 km) (Schmidt and Poli, 1998). In the oceanic crust, water released is associated with mineral dehydration such as amphibole, chlorite, zoesite, and lawsonite. After amphibole destabilization (>2.2–2.4 GPa, i.e. about 70 km depth), the water content in the oceanic crust decreases sharply, falling below 1 wt%. Above 3–4 GPa (i.e. about 100 km depth), lawsonite is the last hydrous mineral to remain stable (Schmidt and Poli, 1998). It is noteworthy that lawsonite exhibits an impressive water storage capacity of 11.3 wt%. Hydrated minerals in the lithospheric mantle include serpentine (antigorite, lizardite, and chrysotile) and chlorite at greater depths. Lizardite and chrysotile form at shallow depths, remaining stable below 300 °C. Deeper, we find antigorite, which is stable between 320 °C and 650 °C (i.e. about 90 km depth) (Schwartz et al., 2001, 2013). Finally, water can be transported even deeper through chlorite, which remains stable up to 800 °C (Van Keken et al., 2011).

Water released from the slab can hydrate the mantle wedge, by forming minerals such as serpentine, chlorite, and talc. If
<table><caption>**Table 1.** Rheological parameters.</caption><thead><tr><th></th><th>Sediment</th><th>Oceanic crust</th><th>Continental crust</th><th>Mantle</th><th>Weak zone</th></tr></thead><tbody><tr><td><strong>Friction</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>C<sub>0</sub>–C<sub>∞</sub> (MPa)</td><td>2–1</td><td>2–1</td><td>20–2</td><td>20–2</td><td>1–1</td></tr><tr><td>φ<sub>0</sub> – φ<sub>∞</sub> (°)</td><td>2–0.1</td><td>2–0.1</td><td>20–10</td><td>20–10</td><td>0.1–0.1</td></tr><tr><td>ε<sub>min</sub> – ε<sub>max</sub></td><td>0.5–2</td><td>0.5–2</td><td>0</td><td>0.5–2</td><td>0.5–2</td></tr><tr><td><strong>Pseudo-plasticity</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>σ<sub>max</sub> (MPa)</td><td>400</td><td>400</td><td>400</td><td>400</td><td>400</td></tr><tr><td><strong>Dislocation</strong></td><td>Wet quartzite [1]</td><td>Diabase [2]</td><td>Dry quartzite [3]</td><td>Dry Olivine [4]</td><td>Dry Olivine [4]</td></tr><tr><td>n</td><td>4.00</td><td>3.05</td><td>2.40</td><td>3.50</td><td>3.50</td></tr><tr><td>A<sub>ds</sub> (Pa<sup>-n</sup>s<sup>-1</sup>)</td><td>1.1×10<sup>-28</sup></td><td>6.31×10<sup>-20</sup></td><td>2.5×10<sup>-20</sup></td><td>2.42×10<sup>-16</sup></td><td>2.42×10<sup>-16</sup></td></tr><tr><td>Q<sub>ds</sub> (J mol<sup>-1</sup>)</td><td>2.23×10<sup>5</sup></td><td>2.76×10<sup>5</sup></td><td>1.56×10<sup>5</sup></td><td>5.40×10<sup>5</sup></td><td>5.40×10<sup>5</sup></td></tr><tr><td>V<sub>ds</sub> (m<sup>-3</sup> mol<sup>-1</sup>)</td><td>10×10<sup>-6</sup></td><td>10×10<sup>-6</sup></td><td>10×10<sup>-6</sup></td><td>14×10<sup>-6</sup></td><td>10×10<sup>-6</sup></td></tr><tr><td><strong>Diffusion</strong></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>A<sub>df</sub> (Pa s<sup>-1</sup>)</td><td></td><td></td><td></td><td>1.92×10<sup>-11</sup></td><td></td></tr><tr><td>Q<sub>df</sub> (J mol<sup>-1</sup>)</td><td></td><td></td><td></td><td>3.00×10<sup>5</sup></td><td></td></tr><tr><td>V<sub>df</sub> (m<sup>3</sup> mol<sup>-1</sup>)</td><td></td><td></td><td></td><td>4.5×10<sup>-6</sup></td><td></td></tr></tbody></table>

Gleason and Tullis (1995) [1]; Burov and Guillou-Frottier (1999) [2]; Ranalli and Murphy (1987) [3]; van Hunen et al. (2000) [4].

the water is released at greater depths, it will trigger partial melting of the mantle (Hyndman and Peacock, 2003). Several studies (Van Keken et al., 2011; Hacker, 2008; Gies et al., 2024) have predicted the dehydration of subducting slabs across various subduction zones worldwide using thermodynamic modeling combined with a kinematic-dynamic approach. Their results reveal a strong dependence on the geothermal gradient, with significantly different dehydration timings in warm versus cold subduction zones.

In this study we investigate how slab dip angle influences the water budget, specifically evaluating the contributions of sediment, crust, and lithospheric mantle to melting and hydration of the mantle wedge. Observations related to volcanism in the Pampean flat slab segment suggest that flat subduction is a dynamic process evolving over time. Consequently, a purely thermodynamic approach appears insufficient for studying the water budget in this context. To address this, we have implemented water circulation within the pTatin2D thermomechanical code (May et al., 2014, 2015, Jourdon et al., 2018; Larvet et al., 2022) and conducted simulations for both steep and flat subduction scenarios. Finally, we discuss the implications of the water budget for volcanism and H₂ production in the mantle wedge.

## 2 Method

The pTatin2D code is a thermomechanical modelling software specifically designed for studying long-term geo-dynamic processes. It has been previously utilized and validated for analysing subduction zones (Larvet et al., 2022, 2023). Our study presented an opportunity to incorporate fluid circulation and simulate the partial melting of mantle and oceanic crust. In the following section, we outline the key equations governing pTatin2D, provide an explanation of how water circulation was implemented, and describe our simulations.

### 2.1 Governing equations

The pTatin2D employs the finite element method on an arbitrary Lagrangian-Eulerian grid combined with the material point method. At each time step, it solves the Stokes equations for velocity (v) and pressure (P) for static equilibrium (1), assuming incompressible fluid (Eq. (2)) for velocity. The velocity (computed on the grid nodes) is then projected onto material points which are advected to carry information about rocks and especially rheology (Tab. 1).

$$ \nabla \cdot \eta (\nabla v + (\nabla v)^t) - \nabla P + \rho g = 0 \quad (1) $$

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

In equation (1), the coefficients include the effective viscosity (η), gravity acceleration vector (g), and material density (ρ).

Density is tabulated for each lithology as a function of pressure and temperature conditions (Larvet et al., 2022).

In equation (1), the viscosity (η) is first determined by taking the harmonic average (3) of dislocation creep (4) and diffusion creep (5).

$$ \frac{1}{\eta_v} = \frac{1}{2} \left( \frac{1}{\eta_{ds}} + \frac{1}{\eta_{df}} \right) \quad (3) $$

$$ \eta_{ds} = A_{ds} \frac{1}{\dot{e}^t} \frac{1}{n} - 1e^{-\frac{Q_{ds} + PV_{ds}}{nRT}} \quad (4) $$

$$ \eta_{df} = A_{df}^{-1} e^{-\frac{Q_{df} + PV_{df}}{RT}} \quad (5) $$
A represents the pre-exponential factor, n is the stress exponent, Q stands for the activation energy, V represents the activation volume (as specified in Tab. 1), and the second invariant of the strain rate tensor is denoted as $\dot{e}^{II}$. The result of equation (3) is used to calculate a viscous stress predictor, if the stress exceeds either the limit determined by the Drucker-Prager yield criterion

$$ \sigma_{DP} = C \cos \phi + P \sin \phi, \quad (6) $$

which depends on cohesion (C) and friction angle (φ) or the maximum plastic strength $\sigma_{max}$ fixed at 400 MPa (Watremez et al., 2013), the viscosity is adjusted accordingly to

$$ \eta = \min \left( \eta_v, \frac{\sigma_{\max}}{2\dot{e}^{II}}, \frac{\sigma_{DP}}{2\dot{e}^{II}} \right) \quad (7) $$

The temperature is obtained solving the heat conservation

$$ \frac{\partial T}{\partial t} = -\nabla \cdot (-\kappa \nabla T + \nu \nabla T) + \frac{H}{\rho c_p} \quad (8) $$

for temperature (T), using the thermal diffusivity (κ), the heat capacity (cP), and the heat production (H) computed on markers and projected on the Gauss points of a Q1 FEM grid aligned to the Q2 mechanical grid.

## 2.2 Implementation water and melt

In this study, we consider that (1) the water can be released or stored in the hydrated mineral or be free if the system is saturated; (2) free water movement is governed by a percolation constant and the velocity of the surrounding rocks, (3) melt (and the water it contains) is extracted instantaneously, removing some water from the system and (4) the time to exchange water is smaller than our timestep at the scale of one element (Fig. 2).

### 2.2.1 Water(s)

We refer to different quantities of water. Bulk water is the total amount of water present in the system. It is split into bounded water that is stored in hydrated minerals and free water that is in excess. Note that the mineralogical composition of rocks (sediments, basalt and peridotite) varies depending on the amount of bulk water present. To evaluate bounded water, we utilize precomputed tables (Fig. 3) generated with Perple_X (Connolly, 2005). Perple_X is a software that minimizes Gibbs free energy to predict the stable mineralogical assemblage as a function of thermodynamic variables. Here, we vary P-T (Pressure-Temperature) conditions for a set of fixed bulk rock chemical composition compiled by Li et al. (2019) and listed in Table 2 to obtain first mineralogical assemblage and therefore their hydrated mineral composition. At saturation, the stable mineralogical assemblage exhibits the largest capacity to store bounded water. We define water capacity as the bounded water computed at saturation. The bulk water in excess of water capacity is free to be transported.

### 2.2.2 Free water transport

For water transport, we opted for the kinematic approach developed in Gerya and Meilick (2011) to model oceanic slab dehydration and mantle hydration. As a result, water released by the subducting plate is stored on free water markers that ascend to the surface due to buoyancy at a velocity that is set equal to the velocity of the surrounding rocks plus a user-defined constant vertical percolation velocity arbitrarily set to 5 cm yr⁻¹. As we assume that water within the mantle does not rely on preferential paths like faults to facilitate its movement, water is distributed uniformly throughout the mantle, without exhibiting localized high concentrations.

### 2.2.3 Interactions with partial melting

Partial melting of the hydrated mantle and hydrated oceanic crust are also taken into consideration. For the mantle, the capacity to melt of markers is computed using an analytical solution based on laboratory experiments (Katz et al., 2003). To estimate the water capacity of peridotites, we combine the water capacity associated with hydrated minerals (low temperatures) predicted by Perple_X with the ability of the magma to retain water (high temperatures), as predicted by Katz et al. (2003) (Figs. 3d, 3e, 3f). We also forecast the percentage of melting (Figs. 3h, 3i). Initially, mantle peridotites are regarded as lherzolites. Following Larvet et al. (2022), if the melting percentage of a lherzolite surpasses 22%, the rock is deemed harzburgite, which can no longer undergo melting. As the melting of the mantle is influenced not only by P-T conditions but also by the quantity of available water (free and bounded), we therefore refer to various tables derived from different initial bulk water values (1 wt%, 2 wt%, 12 wt%) taking the one that is the closest to the amount of water carried by the element. For the melting of the hydrated oceanic crust, the amount of melting increases linearly (from 0% to 100%) between the hydrated solidus (Poli and Schmidt, 2002) and hydrated liquidus (Green, 1982) of the basalt. Finally, at each time step, the melt is extracted by resetting the water content within the melted rocks to zero. This way, some water is removed from the system.

### 2.2.4 Water conservation at element level

For each element, we use equation (9) to predict if it can store more water ($M_{H_2O} < 0$) or if it has exceeded its maximum capacity ($M_{H_2O} > 0$).

$$ M_{H_2O} = \frac{S}{n} \sum_{k=1}^{n} \rho_k (bw_k - wc_k) \quad (9) $$

where n represents the number of markers in the element, S denotes the surface area of the element [m²], and $\rho_k$, $bw_k$ and $wc_k$ respectively represent the density [kg.m⁻³], the bounded water value [wt%] and the water capacity [wt%] of a marker.

Even if an element is able to store more water ($M_{H_2O} < 0$), some markers inside might already be full ($bw > wk$). In such
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Explanatory diagram of water transport at element level.</figcaption></figure>

cases, the excess water that these markers cannot hold is redistributed within the element to markers that possess the capacity for further water storage ($bw < wk$). If an element has surpassed its ability to store water because of the increasing $P-T$ conditions ($M_{H_2O} > 0$), it undergoes dehydration. At this stage, the bounded water values for the markers associated with that element are adjusted to match their respective water capacities. Additionally, we introduce a new marker which carries information about the mass of free water. These markers move within the mantle according to Section 2.2.2 until they encounter an element capable of storing water ($M_{H_2O} > 0$).

### 2.2.5 Monitoring hydration and melting

To evaluate the intensity of hydration in the mantle wedge in our model, we sum the amount of water absorbed by the markers in the asthenospheric mantle and continental lithosphere every 100-time steps (approximately 1 Ma) under hydration conditions. This value is then normalized over 1 Ma. The intensity of melting is assessed in the same way, but this time we consider the markers of the asthenospheric mantle, the continental lithosphere, and also the oceanic crust (which can

also melt) under melting conditions. The values thus obtained are shown on the graphs in Figure 4.

## 2.3 Model set up

Statistical analysis conducted on subduction zones worldwide reveals a correlation between the velocity of the overriding plate and the dip of subduction (Jarrard, 1986; Lallemand et al., 2005). This correlation has been reinforced by numerical studies (van Hunen et al., 2004; Heuret et al., 2007; Huangfu et al., 2016), showing that high absolute velocities of the overriding plate promote flat subductions, as well as positive buoyancy of the subducting plate (van Hunen et al., 2004; Espurt et al., 2008), particularly associated with young plates or the presence of oceanic plateaus. In South America, the observed variations in subduction dip over time are likely attributable to subduction of oceanic plateaus (Martinod et al., 2010). Finally, eastward flow of the asthenospheric mantle (Ficini et al., 2017; Gerya, 2022) can also play a role in decreasing subduction dip (Huangfu et al., 2016). Here, the objective is not to investigate the factors leading to flat subduction, but rather to elucidate how changes
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Perple_X tables for water capacity and melt prediction. x-axis represents temperature (°C) and y-axis represents pressure (GPa). Red line represents hydrated solidus. (a) Sediment water capacity; (b) and (c) oceanic crust (basalt compound) water capacity and melt respectively; (d)–(f) mantle (peridotite compound) water capacity with 1 wt%, 2 wt%, and 12 wt% bulk water content respectively; (g)–(i) mantle melt with 1 wt%, 2 wt%, and 12 wt% bulk water content respectively.</figcaption></figure>

<table><caption>**Table 2.** Rock compositions [wt%] used for water capacity calculations.</caption><thead><tr><th></th><th>MgO</th><th>FeO</th><th>CaO</th><th>Al<sub>2</sub>O<sub>3</sub></th><th>Na<sub>2</sub>O</th><th>K<sub>2</sub>O</th><th>SiO<sub>2</sub></th><th>Reference</th></tr></thead><tbody><tr><td>Sediment</td><td>2.59</td><td>5.43</td><td>6.21</td><td>12.42</td><td>2.53</td><td>2.53</td><td>61.09</td><td>Plank and Langmuir (1998)</td></tr><tr><td>Basalt</td><td>7.57</td><td>10.42</td><td>11.37</td><td>14.68</td><td>2.78</td><td>0.0</td><td>50.4</td><td>Gale et al. (2013)</td></tr><tr><td>Peridotite</td><td>37.44</td><td>7.47</td><td>3.18</td><td>4.03</td><td>0.33</td><td>0.0</td><td>45.55</td><td>Hart and Zindler (1986)</td></tr></tbody></table>
<figure><img src="image_5.png" /></figure>

**Fig. 4.** Amount of water utilized in mantle hydration vs. melting in steep subduction simulations (left, M1 and M5) and flat-slab simulations (right, M2, M3, M6). The green curve represents the convergence rate. The red curve illustrates the amount of water consumed during melting processes, while the blue curve shows the amount of water consumed by mantle hydration. The different background colors represent the different stages of flat subduction simulations: Gray area) Initiation of subduction; Blue area) The dip angle is still relatively high and only mantle melting is observed; Yellow area) Initiation of flat subduction and onset of oceanic crust melting; Green area) Flat subduction is established: the slab lies beneath the continental lithosphere.

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

**Figure 5.** Model setup. We varied mantle wind ($v_{wind}$), overriding plate velocity ($v_{ov}$), and sediment thickness, highlighted in red, in order to obtain a range of flat (M2, M3, M4, M6) and steep (M1, M5) subduction scenario. See Table 3. Subducting plate velocity ($v_{sub}$), varies with time: 6 cm yr$^{-1}$ for the first 10 Ma to initiate subduction, then decreases linearly until reaching 3 cm yr$^{-1}$ from 15 Ma onwards.

in the dip of the subducting plate influence hydration and melting processes within the mantle wedge. To achieve this, we conducted simulations where we varied the velocity of the upper plate and the motion of the asthenosphere, in order to obtain two steep subduction simulations (M1, M5) and three flat subduction simulations (M2, M3, M4). Variable parameters are highlighted in red in Figure 5 and listed in Table 3.

The initial and boundary conditions are constrained by the range of values observed in the Nazca plate subduction below South America between the paleo-flat-slab of Puna and the current flat-slab of Pampean Ranges (Ramos and Folguera, 2009; Martinod et al., 2020) (20°S–30°S) over the past 20 Ma.

During that period, the age of the Nazca Plate in the study area has varied between 60 and 30 Ma (Fig. 1). In our simulations, we therefore used an oceanic lithosphere consistent with an age of 45 Ma. The computational domain (Fig. 5) has dimensions of 2600 km in length and 900 km in depth. In our simulations, the negative Clapeyron slope associated with the ringwoodite-to-perovskite phase transition at a depth of 670 km induces positive buoyancy in the subducting slab relative to the surrounding asthenospheric mantle. This buoyant effect hinders the vertical descent of the slab, preventing it from penetrating beyond 670 km and from reaching the base of the model domain (Larvet et al., 2022). The grid is discretized into
<table><caption>**Table 3.** Varying parameters across models.</caption><thead><tr><th></th><th>Name</th><th>Sediment thickness (km)</th><th>Continental plate velocity (cm yr<sup>-1</sup>)</th><th>Mantle wind (cm yr<sup>-1</sup>)</th><th>Initial serpentinization</th></tr></thead><tbody><tr><td>M1</td><td>Steep</td><td>4</td><td>1</td><td>0</td><td>Yes</td></tr><tr><td>M2</td><td>Flat</td><td>4</td><td>3</td><td>5</td><td>Yes</td></tr><tr><td>M3</td><td>Flat (no initial serpentinization)</td><td>4</td><td>3</td><td>5</td><td>No</td></tr><tr><td>M4</td><td>Flat (static continental plate)</td><td>4</td><td>from 3 to 0 <sup>*1</sup></td><td>5</td><td>Yes</td></tr><tr><td>M5</td><td>Steep (500 m sediment thick)</td><td>0.5</td><td>1</td><td>0</td><td>Yes</td></tr><tr><td>M6</td><td>Flat (500 m sediment thick)</td><td>0.5</td><td>3</td><td>5</td><td>Yes</td></tr></tbody></table>

480 by 96 elements and further refined to achieve a maximum resolution at the subduction interface. The domain is divided into two parts: an oceanic section (the Nazca Plate) on the left and a continental section (the South American Plate) on the right. At the oceanic-continental boundary, a zone of very low friction (weak zone) facilitates the initiation of subduction. The initial angle between the two plates (slope angle) is set equal to 20°.

The subduction process between the Nazca Plate and South America is marked by the substantial absolute motion of the oceanic plate, ranging from 5 to 10 cm yr⁻¹ (Sdrolias and Müller, 2006), and the dynamic nature of the continental plate, which is moving towards the oceanic plate at a velocity of approximately 1 to 3 cm yr⁻¹ (Sdrolias and Müller, 2006). Kinematic boundary conditions are depicted in Figure 5. The lateral boundaries of the domain are permeable (material can leave or enter the domain) with imposed horizontal velocities that are designed to conserve volume in the domain. We fix the velocity acting on the oceanic lithosphere and vary the velocity of the overriding plate. We can also add a translation component between the asthenosphere and the surface that we named ‘mantle wind’. We impose free slip condition at the base of the model and a free surface that capture the topography at the top boundary. Thermal boundary conditions assume insulation of the lateral boundaries and fixed temperature at the top and the bottom of the domain.

The oceanic crust consists a sediment layer (4 km or 0.5 km according to the simulation, Tab. 3) overlaying a 7 km thick basalt layer. The oceanic lithosphere geotherm follows:

$$T_{oce} = \frac{10(T_{Lo} - T_{So})}{9} \operatorname{erf}\left(\frac{y}{2\sqrt{\kappa t}}\right), \quad (10)$$

according to a half space cooling model (Turcotte and Schubert, 2002). In equation (10), $T_{So}$ represents the temperature at the surface (0°C), y denotes the depth (in meters), κ represents the thermal diffusivity coefficient in the continental lithosphere (1e-6 m s⁻¹), and t represents the age of the oceanic plate (in seconds). $T_{Lo}$ corresponds to the temperature of the mantle of the Asthenosphere. Using equation (10), Turcotte and Schubert (2002), the depth of the Lithosphere Asthenosphere Boundary (LAB) is

$$y_{Lo} = -2.32\sqrt{\kappa t}. \quad (11)$$

In our simulations, the continental crust is composed of a 30 km granite layer and the LAB is fixed at $y_{LAB}$=100 km. The thermal structure of the continental region is influenced by a radiogenic heat production that is deemed to decrease exponentially with depth (Lachenbruch, 1968). Assuming heat conservation for a conductive lithosphere at steady state, the temperature of the continental lithosphere ($T_{cont}$) is determined by the equation (12).

$$T_{cont} = T_{Sc} + \left( \frac{T_{LAB} - T_{Sc}}{y_{LAB}} \right) y + \frac{h_r^2 S_o}{k} \left( 1 - \exp \left( -\frac{y}{h_r} \right) \right) \quad (12)$$

In the equation, $T_{Sc}$ represents the surface temperature (20°C), and k represents the thermal conductivity (3 J s⁻² m⁻¹ K⁻¹), $h_r$ denotes the characteristic length scale of decrease for heat production (10 km).

The asthenospheric geotherm is characterized by a gradient of 0.5 °C km⁻¹ and is anchored at a specific point determined by the pressure and temperature conditions associated with the destabilization of Ringwoodite and the formation of Bridgmanite (1600 °C at 670 km).

At the onset of the simulations, the quantity of bounded water in the rocks of the oceanic lithosphere is initialized by assigning each marker a value representing the percentage of bounded water [wt%] (i.e. hydroxyl group -OH within the structure of some minerals) before subduction (i.e. hydrothermalized seafloor). The initial amount of bounded water relies on proposal from Faccenda (2014). Sediments have the highest water content (7 wt%). Furthermore, the water content of the oceanic crust decreases with depth, ranging from 3 to 0.5 wt%. In central Chile, as deduced from the interpretation of seismic data (Grevesmeyer et al., 2007; Kopp et al., 2004), the uppermost section of the oceanic lithosphere exhibits signs of water alteration. Therefore, 10 % serpentinization (i.e. approximately 1.2 wt%) in the first 5 km, is considered (Faccenda, 2014).

# 3 Results

## 3.1 Steep vs. flat subduction

Figure 6 illustrates the evolution of water circulation in the case of steep subduction (M1) compared to flat subduction (M2). In both simulations, water is released from the slab through the successive destabilization of hydrated minerals. In

<footnote>*1 3 cm yr⁻¹ for the first 10 Ma, then decreases linearly until reaching 0 cm yr⁻¹ from 15 Ma onwards.</footnote>
<figure><img src="image_7.png" /></figure>

**Fig. 6.** Steep vs. flat subduction over time. The boxed value is the amount of water consumed in mantle hydration during the last time step (normalized to 1 Ma). For steep subduction, the hydration zone (in green) is limited to a small area before the volcanic arc (triangle). The water released by the subducting slab, which is not consumed in mantle hydration, is consumed in mantle melting (in red). On the other hand, in flat subduction, the hydration zone is much more extensive. Additionally, it increases over time, as the volcanic arc migrates landward. Oceanic melting (orange) is also observed.

the case of steep subduction, water released at shallow depths is involved in the hydration of the upper plate lithospheric mantle (depicted in green). At higher depths, water is consumed by partial melting (depicted in red) of the asthenosphere mantle (Fig. 6a-c). This distribution remains stationary over time, on the contrary, it evolves rapidly in the flat slab simulation. During the first stage of the simulation, the slab dip is still relatively high, and, as in the steep simulation, the water released by the slab hydrates the upper plate lithosphere or participates in the melting of the asthenosphere mantle wedge. From 13 Ma (Fig. 6d), the slab dip angle reduces, and the slab flattens under the upper plate lithosphere. This change in subduction dynamics generates a shift in the temperature conditions of the slab at depth. While the quasi-absence of an asthenosphere mantle wedge promotes mantle wedge hydration (at the expense of mantle melting), the relatively high temperature in the slab leads to the melting of the oceanic crust (depicted in orange, Fig. 6d-f). After 15 Ma of flat slab, the subduction zone cools down, the partial melting of the hydrated crust decreases, and mainly the hydration process remains (Fig. 6f).

## 3.2 Spatial distribution of hydration/melting processes

The green area shows the hydrated portion of the mantle since the beginning of the simulation, and the purple area emphasizes the location of active production of hydrated mantle (i.e. mantle hydrated during last time step, Fig. 6). It clearly appears that the landward boundary of the active hydration area corresponds to the onset of magmatic activity.

In steep subduction (M1), the hydration of mantle rock is restricted to a zone of 30 km within the forearc portion of the upper plate lithosphere (20 km after 13 Ma, 28 km after 19 Ma, and 29 km after 25 Ma). In the flat subduction simulation (M2), the lateral extent of hydrated mantle increases over time (97.7 km after 13 Ma, 110 km after 19 Ma, and 298 km after 25 Ma) due to landward migration of the magmatic arc over time. The coloured bars (red for mantle melting and orange for oceanic crust melting) are the vertical projection of the melting zone on the surface, delineating areas where volcanism may potentially be active. In the M2 simulation, after 19 Ma (Fig. 6e), flat subduction is established, i.e. the slab is located just beneath the continental lithosphere. From there, two distinct oceanic crust melting zone, triggered by different processes, are observed (e.g. Fig. 7b, Zone A and Zone B). The geotherm in flat subductions is warmer than in steep subductions, therefore oceanic crust hits the melting conditions (defined by the hydrated solidus of basalt) before it completely loses its bounded water resulting in oceanic crust melting. At the same time, oceanic crust melting is fuelled by the final dehydration of sediments; Zone B) Chlorite breakdown hydrates the oceanic crust above it, lowering melting conditions and triggering oceanic crust melting. Note that the two melting zones are not always distinguishable in all our figures depicting flat subductions. After 25 Ma, only zone A is active for M4 (Fig. 7c), while only zone B is active for M6 (Fig. 8b). This variation in activity of the melting zones does not seem to be linked to subduction parameters but appears to be associated with temporal variations related to punctual numerical artefact. For instance, in M6, both melting zones are active at 23 Ma (see supplementary material videos).
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** (a) No initial serpentinization of the slab: reduction in mantle hydration and disappearance of the second melting zone. (b) Zones A and B both show oceanic plate melting, but triggered by different processes (see Sect. 3.2). Zones I, II, III refer to Section 4.2.2. (c) The extent of hydration zone depends on the velocity of the overriding plate. See the legend Figure 6.</figcaption></figure>

To further investigate the migration of volcanic arcs and the lateral extent of mantle hydration zone during flat subduction, we carried out simulation M4, which differs from M2 in the velocity imposed on the overriding plate (Tab. 3). In simulation M4, we applied a large velocity (3 cm yr⁻¹) to the continental plate for 10 Ma to initiate the reduction of the slab dip angle. Then, the velocity of the continental plate gradually decreases, reaching zero from 15 Ma. As the convergence velocity is reduced, flat slab segment is shorter for M4 than for M2 (Fig. 7b, c). In both simulations, slab dehydration starts at a similar distance from the accretionary prism: 261 km for M4 and 285 km for M2. However, after 25 Ma, the size of the hydrated zone is almost three time larger for M2 (298 km, Fig. 7b) than for M4 (111 km, Fig. 7c).

To monitor the evolution of arc migration we measure the distance between the prism and beginning of the arc (D<sub>PA</sub>) (Figs. 6,7,9). In the steep simulation, this distance is short and remains constant through time (about 250 km). On the contrary, in the flat slab simulation with constant upper plate v plate velocity (M2), it keeps growing over time: 322 km after 13 Ma (Fig. 5d), 399 km after 19 Ma (Fig. 6e), and 583 km after 25 Ma (Fig. 5f), while the distance is reduced (499 km for M2 versus only 365 km for M4, Figure 7b, c) in the simulation where the upper plate stops moving after 15 Ma (M4). This implies that arc migration and the extent of the hydrated lithosphere portion of the upper plate depend on the length of the flat slab segment which, in turn, depends on the overriding plate velocity.

## 3.3 Water budget

In our simulations, water can (1) be released in the prism at very shallow depths, (2) hydrate the mantle, (3) triggers melting. Figure 9 provides a comparison between the amount of water utilized in mantle hydration and melting for steep (Fig. 9a) and flat (Fig. 9b) subduction simulations over time. Convergence velocity is also represented.

For steep subduction scenarios (M1), the melting rate increase during the first 8 Ma and then decreases from 15 Ma in response to the reduced water flux resulting from the deceleration of the slab (Fig. 9a). The rate of water hydrating the mantle remains relatively stable, around 2 Mt Ma⁻¹ m⁻¹. Conversely, for M2, a shift in geothermal conditions occurs from the reduction of the slab dip angle initiation at 13 Ma. This results in increased water release by the slab before reaching the hydration/fusion limit, and less afterwards. Consequently, the flat subduction scenario (M2) shows higher water consumption in mantle hydration reactions compared to steep subduction (M1), which amplifies over time (10.1 Mt Ma⁻¹ m⁻¹ at 13 Ma, 9.5 Mt Ma⁻¹ m⁻¹ at 19 Ma, 20.1 Mt Ma⁻¹ m⁻¹ at 25 Ma). Similarly, in flat subduction context, the decrease in melting intensity is more pronounced (Fig. 9). Indeed, at 25 Ma, melting intensity is higher for M1 than for M2 (14.8 Mt Ma⁻¹ m⁻¹ vs. 9.7 Mt Ma⁻¹ m⁻¹), despite the water flux is more important for M2.

To study the influence of chlorite breakdown, we run simulation M3 (Fig. 7a). This simulation is similar to M2, but without considering the initial water alteration of the oceanic
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Flat vs. steep subduction with 500 m sediment thickness. The reduction in the quantity of sediment subducting strongly reduces mantle wedge hydration.</figcaption></figure>

lithosphere (Tab. 3). Therefore, as expected, the second melting zone (Zone B), which is triggered by chlorite breakdown, is not observed for M3. This observation is persistent over time and cannot be associated with a numerical artefact. Despite a higher amount of water entering subduction when the lower plate lithosphere is serpentinized (M2) comparing to M3, the quantity of water reacting with the mantle after 25 Ma simulation is similar in simulations M3 (22.7 Mt Ma⁻¹ m⁻¹) and M2 (20.1 Mt Ma⁻¹ m⁻¹) which show that the water release by chlorite breakdown does not participate to mantle hydration process but only to melting. In the case of steep subduction, the water goes into the asthenosphere mantle and generate the melting. In the case of flat subduction, water goes into the anhydrous portion of the basaltic crust above it and lead to oceanic crust melting. Hence, lower plate serpentinization does not influence mantle hydration. When flat subduction is established, the intensity of melting is much greater when the lithosphere is initially serpentinized (Fig. 9b): after 25 Ma of subduction, the amount of water used in melting is only 2.2 Mt Ma⁻¹ m⁻¹ for M3 versus 9.7 Mt Ma⁻¹ m⁻¹ for M2. This shows that the amount of magma generated in zone A is low, while it is higher in zone B.

During subduction, sediments from the oceanic crust can either be accreted into the prism or be dragged deeper into the slab, where they participate in hydration/melting processes. To highlight their role in the water budget we conducted simulations M5 (steep subduction) and M6 (flat subduction), which are respectively the same as M1 and M2 but with only 500 m instead of 4 km of sediment initially deposited on the oceanic crust (Fig. 8). For M1 and M2, a significant proportion of the sediments are subducted, whereas for M5 and M6, they are exclusively accreted into the prism. For M5 (steep subduction) at 25 Ma, the surface area of the hydrated mantle wedge (green area) is very small. Similarly, active hydration (purple area) is almost non-existent: only 0.4 Mt Ma⁻¹ m⁻¹ of water hydrates the mantle, compared to 1.8 Mt Ma⁻¹ m⁻¹ for M1. This means that in our simulations under steep subduction conditions, it is almost exclusively the sediments that hydrate

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Horizontal distance between the eastern boundary of the accretionary prism and the volcanic arc ($D_{PA}$). Flat subduction scenarios show a landward migration of the volcanic arc, while in steep subduction scenario the distance between the prism and the arc is constant over time.</figcaption></figure>

the mantle wedge. In the context of flat subduction, mantle hydration is much more significant (Fig. 4), indicating that changes in geothermal conditions allow the crust to hydrate the mantle. However, mantle hydration is still considerably less intense than for M2 (6.0 Mt Ma⁻¹ m⁻¹ vs. 20.1 Mt Ma⁻¹ m⁻¹). This shows that even in a flat subduction context, sediments play an important role in the hydration budget of the mantle wedge.

# 4 Discussion

## 4.1 Steep subduction

Our results for steep subduction (M1 and M5) are consistent with thermodynamic studies (Van Keken *et al.*, 2011; Hacker, 2008; Gies *et al.*, 2024). The sediments
primarily dehydrate at shallow depths, contributing to mantle
hydration, while the oceanic lithosphere dehydrates at greater
depths, driving partial mantle melting (Van Keken et al., 2011).
However, our study has highlighted the crucial role of
sediments in the hydration of the mantle wedge. For steep
subduction simulations, the mantle wedge is almost exclu-
sively hydrated by the sediments. These results should be
interpreted with caution as they can vary significantly
depending on the rock composition considered (Hacker,
2008) or the initial amount of water present in the subducting
slab (Van Keken et al., 2011). These parameters vary greatly
depending on the specific case studied. Here, we have used
global subducting sediment (Plank and Langmuir, 1998).

The results indicate that, in the case of steep subductions (M1), the mantle hydration zone is located between the prism and the magmatic arc (see Fig. 6). However, this zone is limited to 30 km before the arc and does not form a continuous connection with the prism. Similar studies (Arcay et al., 2005; Gerya and Meilick, 2011; Menant et al., 2020), using the same implementation method (kinematic approach) for water transport as in this study (see Gerya, 2022 for review) give results in agreement with ours. For example, Menant et al. (2020) showed a mantle hydration zone limited to 40 km with the magmatic arc with a D<sub>PA</sub> around 200 km. They also highlighted the uncertainties associated with the percolation constant. While this constant is geologically poorly con- strained, it exerts an influence on the results. A lower value accentuates the relative significance of the surrounding rocks in determining the velocity of free water within the mantle. Within the mantle, as rocks migrate towards the subduction trench (due to the trenchward motion of continental lithosphere), a diminished percolation constant corresponds to a larger extension of the mantle hydration zone towards the prism. Simulations from Gerya and Meilick (2011) shows a continuous hydration of the mantle between the prism and the volcanic arc. This distribution contrast comes from the initial amount of water present in the oceanic plate, which is substantially greater than in our initial condition constrained by Faccenda (2014). During subduction, as P–T conditions increase in the slab, the ability to retain water diminishes. In the case of Gerya and Meilick (2011), as the slab contains more water, a part of the water released by dehydration below the prism goes into the lithospheric mantle of the upper plate, whereas in our simulations, it is totally recycled within the slab to form new hydrated minerals stable at higher P–T conditions.

More physics-based implementations have also been
explored (see Gerya, 2022 for review). Angiboust et al.
(2012), utilized the Darcy equation to promote fluid migration
within deformation zones characterized by increased porosity.
This method highlights specific paths for water migration,
leading to localized hydration—unlike our simulations, where
hydration is homogeneously distributed. To account for
porosity variations, several authors have also adopted the
two-phase flow approach (Cerpa et al., 2019; Wang et al.,
2019; Wilson et al., 2014). Although this method doesn't
predict dehydration and hydration processes, it illustrates the
migration of free water into the mantle wedge. The water
released by the slab ascend along preferential path and
especially along the subduction interface. It results in
distributed hydration between the prism and the volcanic

arc. Geological and geophysical studies conducted in many
steep subductions around the world (Bebout and Barton, 1989;
Hyndman and Peacock, 2003; Myers et al., 1998) have also
demonstrate this distributed hydration. This implies that H₂
production in the mantle wedge would extends from the prism
to the arc in a steep subduction context, resulting in more
widespread and therefore less concentrated H₂ production than
in the M1 simulation.

4.2 Flat subduction

Simulations with flat subduction scenarios have shown that the hydration zone is much larger than in steep subductions, and that it increases with time. The extent of hydration zone seems to depend on the velocity of the overriding plate. There are no other studies presenting thermomechanical simulations of mantle hydration in flat subduction contexts but we can compare our results with geological and geophysical observations in the Pampean flat slab.

4.2.1 Magmatism in the Pampean flat slab

Pampean flat slab (Fig. 1), began during early Miocene,
around 18 Ma ago with the reduction of the subduction angle of
the Nazca Plate beneath South America (Litvak et al., 2007).
As observed in our simulations (Fig. 9b), this reduction in slab
dip angle is marked by an landward migration of the volcanic
arc (Kay and Abbruzzi, 1996; Litvak et al., 2007; Ramos et al.,
2002). Ramos and Folguera (2009) demonstrate that volcanism
extended up to 600 km from the prism, in line with our
simulation M2 in Figure 6, where melting is between 400 km
and 800 km approximately.

Initiation of the Pampean flat slab was accompanied by a broadening of the volcanic arc (Jones et al., 2016), which can be seen in the M2 simulation (Fig. 6d) and a decrease in the amount of lava expelled by volcanoes (Ramos and Folguera, 2009) until cessation 4.5 Ma ago (Poma et al., 2017). According to our simulations, the initiation of the flat slab is indeed associated with a decrease in magma production (Fig. 9b). This decrease in magma production, combined with the expansion of the melting zone, results in a drastic reduction in magma concentration. These observations could potentially explain the cessation of active volcanism in the region.

Additionally, from the Pampean flat slab, change in magma
chemistry (increase in the La/Yb) and presence of adakitic-
type volcanoes (Kay and Mpodozis, 2002; Poma et al., 2017,
2023) is observed. Adakitic volcanism can be triggered by
several processes (R.W. Kay and Kay, 2002). While Litvak
et al. (2007), and Poma et al. (2017, 2023) proposed that the
adakitic volcanism of the Pampean flat slab results from
magma contamination by thickened continental crust, many
others generally attribute the presence of adakitic volcanism to
oceanic crust melting at slab tear or slab break-off (Pe-Piper
and Piper, 1994; Omrani et al., 2008; Wu et al., 2022) in
different geodynamic context. Here, like Gutscher et al.
(2000), we propose here that the adakitic signature of the
volcanoes during Pampean flat slab is related to oceanic crust
melting due to flat subduction context. This scenario results in
adakitic signature, without invoking slab tear or slab break-off,
that is consistent with geological observations, without
excluding the possibility of a slight crustal contamination.
## 4.2.2 Mantle hydration

The hydration of the mantle wedge of the Pampean flat slab has been studied using seismic analysis of $V_p$ (P-wave velocity) and $V_s$ (S-wave velocity) (Linkimer et al., 2020; Marot et al., 2014; Wagner et al., 2008). These studies identified three successive zones going from west to east. Zone I is the westernmost zone and begins from the prism. It consists of a hydrated mantle with minerals such as chlorite and serpentine. Next, Zone II is dry mantle. Finally, Zone III, the easternmost zone, is characterized by a sharp reduction in $V_p$ and $V_s$, and an increase in $V_p/V_s$ ratio. Linkimer et al. (2020) identified it as a highly hydrated zone (5 wt%), but $V_p - V_s$ values also fit with melt.

These zones are also present for M2 at 25 Ma (Fig. 7b). Zone I extends up to 583 km from the prism, compared with around 300 km for Linkimer et al. (2020). The extent of this zone depends notably on the velocity of the overriding plate (measuring only 429 km for M4) but also potentially on other parameters that we have not studied, such as the age of the oceanic plate, for example. As with normal subductions, geophysical studies predict hydration up to the prism, which is not the case in our simulation (see Sect. 4.1). Then, Zone II measures 80 km for M4 versus 50 km for Linkimer et al. (2020). This zone of dry mantle is determined by the timing of the hydrated minerals destabilization of the slab. It is located beyond complete dehydration of the oceanic crust and before chlorite breakdown. Zone III is also present in M2 and corresponds to Zone B. In our simulation, it corresponds to melting due to chlorite dehydration, which is consistent with the measured seismic velocity. This means that even in a context of flat subduction with extinct volcanism, magma production can still be active and can even be observed with geophysical measurements. However, the first melting zone (Zone A) is not described by Linkimer et al. (2020). Our models have shown that Zone A produces only a small quantity of magma, which may explain why it is not detected by geophysical measurements.

## 4.3 Implication for deep H₂ generation

The combustion of Dihydrogen ($H_2$) produces energy and water. As such, it is seen as a promising solution to reduce carbon emissions in industry and transportation. Artificial production methods, including electrolysis of water or deriving it from hydrocarbons have been explored, but emit $CO_2$ and/or consume a lot of energy. In contrast, natural sources of $H_2$ (also called white dihydrogen) could offer a more eco-friendly and cost-effective production process (Gaucher et al., 2023; Moretti and Webber, 2021). The natural production of $H_2$ occurs through several geochemical processes including reduction of water via rock alteration. Among water reduction processes, serpentinization is the most widely recognized and is known to generate $H_2$. Serpentinization occurs when olivine minerals hydrate under favourable P-T conditions, resulting in the formation of serpentine, brucite, magnetite, and $H_2$ (McCollom et al., 2022). Antigorite, lizardite and chrystotile are the main serpentine minerals. Under the P-T conditions of the mantle wedge, only antigorite can be stable. While $H_2$ production is well known in the formation of low-grade serpentine (lizardite and chrysolite), $H_2$ production during

antigorite formation is open to debate. Based on a thermodynamic study, Evans (2010) concluded that under HP-HT conditions, iron oxidation and hence $H_2$ production are unlikely. However, in recent studies Boutier et al. (2021) and Vitale Brovarone et al. (2020), have observed fluid inclusions of $H_2$ in exhumed serpentinized mantle rocks from the lower plate of the Alpine paleo-subduction zone. Based on textural relations, the authors propose that the $H_2$ was produced in the antigorite stability field which implies deep production of $H_2$ and suggests that $H_2$ can be produced at mantle wedge P-T conditions. Other hydrated minerals could be present in the mantle wedge, such as chlorite, for example (Manthilake et al., 2016). The $H_2$ production associated with the formation of these minerals has not yet been studied at such P-T conditions. Nevertheless, experimental work conducted by Rough (2011) at 420 °C and 500 bar showed that the formation of chlorite, smectite and talc linked to olivine hydration produces $H_2$. Therefore, it is likely that, just like in the case of serpentinization, these reactions, which are also linked to water reduction, generate $H_2$ through the oxidation of $Fe^{2+}$ present in peridotites under HP-HT conditions.

In summary, recent studies seem to show that mantle wedge hydration reactions could generate $H_2$. There are currently no equations to calculate the amount of generated $H_2$ based on the amount of consumed water in the mantle wedge. However, we can assume that the more water is consumed, the more $H_2$ is produced and use the amount of consumed water as a proxy for $H_2$ production. In a different geodynamic context, Zwaan et al. (2023) and Liu et al. (2023) adopted a similar approach to discuss hydrogen fluxes in different geodynamic contexts, respectively rift inversion in orogens and mid-oceanic ridges. They also aimed to evaluate $H_2$ production, not by explicitly modelling the chemical reactions, but by identifying mantle material that could be hydrated under favourable P-T conditions. Note that their hypotheses concern $H_2$ production during low-grade serpentine formation which is a better thermodynamically constrained redox reaction.

Our results highlight the competition between $H_2$ production (mantle hydration) and melting. In steep subduction scenarios, a higher proportion of water is used in melting compared to flat subductions. This leads to an increased mantle hydration, consequently resulting in greater $H_2$ production in flat subduction. Additionally, in a flat subduction context, the quantity of $H_2$ generated exhibits a progressive increase over time, making flat subduction in an advanced state (such as the 18 Ma Pampean flat) particularly appealing. Lastly, flat subductions with a slow convergence rate present a notable advantage for $H_2$ exploitation, as the limitation of hydration zone extent increases $H_2$ concentration. Nevertheless, note that this study focused on $H_2$ production at depth, but cannot predict surface releases, as we do not model the $H_2$ migration from the deep source to the surface.

Our simulations have highlighted the importance of sediments in mantle hydration and consequently in $H_2$ production. In the context of steep subduction, our results show that mantle hydration is almost exclusively linked to water supplied by the sediments. However, it would be incorrect to conclude that the oceanic crust cannot hydrate the mantle wedge in such a context based on these results. For instance, by considering an initially more hydrated oceanic crust, like in Gerya and Meilick (2011), earlier dehydration
would be observed, allowing the oceanic crust to then hydrate the mantle wedge.

In South America, based on the key parameters identified by our study, the deep H₂ generation potential of the Pampean and Peruvian flat slab segments seem to be similar. On one hand, both subductions are in an advanced state, favouring mantle wedge hydration at the expense of melting: the Pampean flat slab began 18 Ma ago (Litvak et al., 2007), while the Peruvian flat slab started 15 Ma ago (Margirier et al., 2015). On the other hand, the convergence rate is equivalent for both segments: it was 10 cm yr⁻¹ 10 Ma ago and it is now 5 cm yr⁻¹ (Sdrolias and Müller, 2006). However, the amount of subducted sediment appears to differ: by comparing the amount of sediment accreted at the trench and the volume of the prism, (Clift et al., 2009) estimated that 64 Mm³ myr⁻¹ enter subduction at the Peruvian segment, compared to only 27 Mm³ myr⁻¹ for the Pampean segment. This would imply greater H₂ production in Peru. However, in a more recent study, Lallemand et al. (2024) also attempted to estimate the amount of sediment entering subduction by characterizing the thickness of the subduction channel. They found similar values (0.8 km for the Pampean flat slab and 0.9 km for the Peruvian flat slab), which imply a similar amount of subducted sediment and thus a similar deep H₂ generation potential.

We used the amount of water consumed in mantle hydration reactions as a proxy for the released H₂. This approximation allowed us to identify subduction dynamics most favourable for H₂ production. To enhance the precision of our predictions, it seems necessary to adopt a different implementation methodology (such as fluid and melt percolation models) to simulate the fluid migration in the mantle more realistically. This would enable us to better anticipate local water concentrations and, consequently, more accurately predict the phases that form. By gaining a deeper understanding of these phases and knowing the chemical equations governing H₂ production during formation of mantle hydrated minerals, it would be possible to predict the quantity of H₂ generated more precisely. Given the current state of knowledge, the method developed in this study, with its assumed approximations, offers the advantage of providing compelling results by excluding from the interpretation phenomena that are not sufficiently constrained (mantle permeability, equations governing H₂ generation...).

# 5 Conclusions

Our simplified simulations highlight that steep and flat slab scenarios result in different water budgets between melting and mantle hydration. In the first case, water from the slab (sediment, crust and oceanic lithosphere) is mostly released in a warm asthenospheric mantle environment, producing mantle melting. In the second case, water is released in a cold mantle environment (the mantle lithosphere of the upper plate), promoting mantle hydration over a wider and wider area with time and cessation of arc volcanism. Moreover, mantle wedge hydration seems to be particularly sensitive to the amount of subducted sediment, whether in flat or steep subduction. Lower plate serpentinization doesn't affect mantle wedge hydration because the water released by

chlorite breakdown participates only to melting. Finally, flat subduction models predict oceanic crust melting, producing magma with adakitic signature without invoking the occurrence slab break-off nor slab tear.

Using the amount of water consumed by mantle hydration reactions as a proxy for the released H₂, we can state that the hydration of the mantle wedge in the flat subduction case is much larger than in steep subduction. These areas are the most promising for H₂ exploration. Our simulations aim to produce a conceptual model for H₂ release at subduction zones. In order to be more predictive for H₂ exploration, e.g. find surface markers of mantle hydration, dedicated simulations including more geological background in their set up, such as pre-existing weakness zones in the upper plate, would need to be performed in the future.

# Supplementary material

**Supplementary video 1:** Video of M1 simulation. See Figure 6 legend.

**Supplementary video 2:** Video of M2 simulation. See Figure 6 legend.

**Supplementary video 3:** Video of M3 simulation. See Figure 6 legend.

**Supplementary video 4:** Video of M4 simulation. See Figure 6 legend.

**Supplementary video 5:** Video of M5 simulation. See Figure 6 legend.

**Supplementary video 6:** Video of M6 simulation. See Figure 6 legend.

The Supplementary Material is available at https://www.bsgf.fr/10.1051/bsgf/2024026/olm.

# References

Angiboust S, Wolf S, Burov E, Agard P, Yamato P. 2012. Effect of fluid circulation on subduction interface tectonic processes: insights from thermo-mechanical numerical modelling. Earth Planet Sci Lett 357-358: 238-248.

Arcay D, Tric E, Doin M-P. 2005. Numerical simulations of subduction zones. *Phys Earth Planet Inter* 149: 133-153.

Bebout GE, Barton MD. 1989. Fluid flow and metasomatism in a subduction zone hydrothermal system: Catalina Schist terrane, California. *Geology* 17: 976.

Boutier A, Brovarone AV, Martinez I, Sissmann O, Mana S. 2021. High-pressure serpentinization and abiotic methane formation in metaperidotite from the Appalachian subduction, northern Vermont. *Lithos* 396: 106190.

Burov EB, Guillou-Frottier L. 1999. Thermomechanical behavior of large ash flow calderas. *J Geophys Res: Solid Earth* 104: 23081-23109.

Carrillo Ramirez A, Gonzalez Penagos F, Rodriguez G, Moretti I. 2023. Natural H2 emissions in Colombian Ophiolites: first findings. *Geosciences* 13: 358.

Cerpa, N. G., Wada, I., & Wilson, C. R. (2019). Effects of fluid influx, fluid viscosity, and fluid density on fluid migration in the mantle wedge and their implications for hydrous melting. *Geosphere*, 15: 1-23. https://doi.org/10.1130/GES01660.1

Clift PD, Vannucchi P, Morgan JP. 2009. Crustal redistribution, crust-mantle recycling and Phanerozoic evolution of the continental crust. *Earth-Sci Rev* 97: 80-04.

Connolly JA. 2005. Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application
to subduction zone decarbonation. *Earth Planet Sci Lett* 236: 524-41.

Espurt N, Funiciello F, Martinod J, Guillaume B, Regard V, Faccenna C, Brusset S. 2008. Flat subduction dynamics and deformation of the South American plate: insights from analog modeling. *Tectonics* 27: 2007TC002175.

Evans BW. 2010. Lizardite versus antigorite serpentinite: Magnetite, hydrogen, and life (?). *Geology* 38: 879-82.

Faccenda M. 2014. Water in the slab: a trilogy. *Tectonophysics* 614: 1-0.

Ficini E, Dal Zilio L, Doglioni C, Gerya T. 2017. Horizontal mantle flow controls subduction dynamics. *Sci Rep* 7: 7550.

Gale A, Dalton CA, Langmuir CH, Su Y, Schilling J-G. 2013. The mean composition of ocean ridge basalts. *Geochem Geophys Geosyst* 14: 489-18.

Gaucher EC, Moretti I, Pélissier N, Burridge G, Gonthier N. 2023. *The place of natural hydrogen in the energy transition: a position paper*. https://doi.org/10.5281/ZENODO8108239

Gerya T. 2022. Numerical modeling of subduction: state of the art and future directions. *Geosphere* 18: 503-61.

Gerya TV, Meilick F. 2011. Geodynamic regimes of subduction under an active margin: effects of rheological weakening by fluids and melts. *J Metamorp Geol* 29: 7-1.

Gies NB, Konrad-Schmolke M, Hermann J. 2024. Modeling the global water cycle—the effect of Mg-Sursassite and Phase A on deep slab dehydration and the global subduction zone water budget. *Geochem Geophys Geosyst* 25: e2024GC011507.

Gleason GC, Tullis J. 1995. A flow law for dislocation creep of quartz aggregates determined with the molten salt cell. *Tectonophysics* 247: 1-3.

Green, T. (1982). Anatexis of mafic crust and high pressure crystallization of andesite. *American Society of Mechanical Engineers (Paper)*, 465-487.

Grevemeyer I, Ranero CR, Flueh ER, Kläschen D, Bialas J. 2007. Passive and active seismological study of bending-related faulting and mantle serpentinization at the Middle America trench. *Earth Planet Sci Lett* 258: 528-42.

Gutscher M-A., Maury R, Eissen J-P., Bourdon E. 2000. Can slab melting be caused by flat subduction? *Geology* 28: 535-38.

Hacker BR. 2008. H₂O subduction beyond arcs. *Geochem Geophys Geosyst* 9: 2007GC001707.

Hart SR, Zindler A. 1986. In search of a bulk-Earth composition. *Chem Geol* 57: 247-67.

Heuret A, Funiciello F, Faccenna C, Lallemand S. 2007. Plate kinematics, slab shape and back-arc stress: A comparison between laboratory models and current subduction zones. *Earth Planet Sci Lett* 256: 473-83.

Huangfu P, Wang Y, Cawood PA, Li Z-H., Fan W, Gerya TV. 2016. Thermo-mechanical controls of flat subduction: insights from numerical modeling. *Gondwana Res* 40: 170-83.

Hyndman RD, Peacock SM. 2003. Serpentinization of the forearc mantle. *Earth Planet Sci Lett* 212: 417-32.

Jarrard RD. 1986. Relations among subduction parameters. *Rev Geophys* 24: 217-84.

Jones RE, Kirstein LA, Kasemann SA, Litvak VD, Poma S, Alonso RN, Hinton R. 2016. The role of changing geodynamics in the progressive contamination of Late Cretaceous to Late Miocene arc magmas in the southern Central Andes. *Lithos* 262: 169-91.

Jourdon A, Le Pourhiet L, Petit C, Rolland Y. 2018. Impact of range-parallel sediment transport on 2D thermo-mechanical models of

mountain belts: application to the Kyrgyz Tien Shan. *Terra Nova* 30: 279-88.

Katz RF, Spiegelman M, Langmuir CH. 2003. A new parameterization of hydrous mantle melting. *Geochem Geophysics, Geosyst* 4.

Kay RW, Kay SM. 2002. Andean adakites: Three ways to make them. *Acta Petrolog Sin* 18: 303-11.

Kay SM, Abbruzzi J. 1996. Magmatic evidence for Neogene lithospheric evolution of the central Andean "flat-slab" between 30 and 32 S. *Tectonophysics* 259: 15-28.

Kay SM, Mpodozis C. 2002. Magmatism as a probe to the Neogene shallowing of the Nazca plate beneath the modern Chilean flat-slab. *J South Am Earth Sci* 15: 39-7.

Kopp H, Flueh ER, Papenberg C, Klaeschen D. 2004. Seismic investigations of the O'Higgins Seamount Group and Juan Fernández Ridge: Aseismic ridge emplacement and lithosphere hydration. *Tectonics* 23.

Lachenbruch AH. 1968. Preliminary geothermal model of the Sierra Nevada. *J Geophys Res* 73: 6977-6989.

Lallemand S, Heuret A, Boutelier D. 2005. On the relationships between slab dip, back-arc stress, upper plate absolute motion, and crustal nature in subduction zones. *Geochem Geophys Geosyst* 6: 2005GC000917.

Lallemand S, Peyret M, Arcay D, Heuret A. 2024. Accretion versus erosion and sediment transfer balance near the subduction interface. *Comptes Rendus. Géoscience* 356: 1-25.

Larvet T, Le Pourhiet L, Agard P. 2022. Cimmerian block detachment from Gondwana: a slab pull origin? *Earth Planet Sci Lett* 596: 117790.

Larvet T, Le Pourhiet L, Pubellier M, Gyomlai T. 2023. Slab pull driven south China Sea opening implies a mesozoic proto South China Sea. *Geophys Res Lett* 50: e2023GL105292.

Li Z-H., Gerya T, Connolly JA. 2019. Variability of subducting slab morphologies in the mantle transition zone: insight from petrological-thermomechanical modeling. *Earth-Sci Rev* 196: 102874.

Linkimer L, Beck S, Zandt G, Alvarado P, Anderson M, Gilbert H, Zhang H. 2020. Lithospheric structure of the Pampean flat slab region from double-difference tomography. *J South Am Earth Sci* 97: 102417.

Litvak VD, Poma S, Kay SM. 2007. Paleogene and Neogene magmatism in the Valle del Cura region: new perspective on the evolution of the Pampean flat slab, San Juan province, Argentina. *J South Am Earth Sci* 24: 117-137.

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

Manthilake G, Bolfan-Casanova N, Novella D, Mookherjee M, Andrault D. 2016. Dehydration of chlorite explains anomalously high electrical conductivity in the mantle wedges. *Sci Adv* 2: e1501631.

Margirier A, Robert X, Audin L, Gautheron C, Bernet M, Hall S, Simon-Labric T. 2015. Slab flattening, magmatism, and surface uplift in the Cordillera Occidental (northern Peru). *Geology* 43: 1031-1034.

Marot M, Monfret T, Gerbault M, Nolet G, Ranalli G, Pardo M. 2014. Flat versus normal subduction zones: A comparison based on 3-D regional traveltime tomography and petrological modelling of central Chile and western Argentina (29°-35°S). *Geophys J Int* 199: 1633-1654.

Martinod J, Gérault M, Husson L, Regard, V. 2020. Widening of the Andes: an interplay between subduction dynamics and crustal wedge tectonics. *Earth-Sci Rev* 204: 103170.
Martinod J, Husson L, Roperch P, Guillaume B, Espurt N. 2010. Horizontal subduction zones, convergence velocity and the building of the Andes. *Earth Planet Sci Lett* 299: 299-309.

May DA, Brown J, Le Pourhiet L. 2014. pTatin3D: High-performance methods for long-term lithospheric dynamics. SC'14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, 274-284.

May DA, Brown J, Le Pourhiet L. 2015. A scalable, matrix-free multigrid preconditioner for finite element discretizations of heterogeneous Stokes flow. *Comput Methods Appl Mech Eng* 290: 496-523.

McCollom TM, Klein F, Ramba M. 2022. Hydrogen generation from serpentinization of iron-rich olivine on Mars, icy moons, and other planetary bodies. *Icarus* 372: 114754.

Menant A, Angiboust S, Gerya T, Lacassin R, Simoes M, Grandin R. 2020. Transient stripping of subducting slabs controls periodic forearc uplift. *Nature Commun* 11: 1823.

Moretti I, Baby P, Alvarez Zapata P, Mendoza RV. 2023. Subduction and hydrogen release: the case of Bolivian Altiplano. *Geosciences* 13: 109.

Moretti I, Webber M. 2021. Natural hydrogen: a geological curiosity or the primary energy source for a low-carbon future. *Renew Matter* 34.

Myers SC, Beck S, Zandt G, Wallace T. 1998. Lithospheric-scale structure across the Bolivian Andes from tomographic images of velocity and attenuation for *P* and *S* waves. *J Geophys Res: Solid Earth* 103: 21233-21252.

Omrani J, Agard P, Whitechurch H, Benoit M, Prouteau G, Jolivet L. 2008. Arc-magmatism and subduction history beneath the Zagros Mountains, Iran: a new report of adakites and geodynamic consequences. *Lithos* 106: 380-398.

Pe-Piper G, Piper DJW. 1994. Miocene magnesian andesites and dacites, Evia, Greece: Adakites associated with subducting slab detachment and extension. *Lithos* 31: 125-140.

Plank T, Langmuir CH. 1998. The chemical composition of subducting sediment and its consequences for the crust and mantle. *Chem Geol* 145: 325-394.

Poli, S., & Schmidt, M. W. (2002). Petrology of subducted slabs. *Annual Review of Earth and Planetary Sciences*, 30: 207-235.

Poma SM, Ramos AM, Litvak VD, Quenardelle SM, Maisonnave EB, Díaz I. 2017. Southern Central Andes Neogene magmatism over the Pampean Flat Slab: implications on crustal and slab melts contribution to magma generation in Precordillera, Western Argentina., Andean Geology 44: 249-274.

Poma S, Ramos A, Litvak VD, Quenardelle S, Maisonnave EB, Rubinstein N. 2023. Orogenic magmatism along the transition Precordillera-Sierras Pampeanas, between 29° and 31° S, Argentina. *J South Am Earth Sci* 123: 104236.

Ramos VA, Cristallini EO, Pérez DJ. 2002. The Pampean flat-slab of the Central Andes. *J South Am Earth Sci*. 15: 59-78.

Ramos VA, Folguera A. 2009. Andean flat-slab subduction through time. *Geolog Soc London Special Publications* 327: 31-54.

Ranalli G, Murphy DC. 1987. Rheological stratification of the lithosphere. *Tectonophysics* 132: 281-295.

Rough MEM. 2011. H2 and chlorite production from an olivine-rich Gabbroic Rock assemblage: A modeling and experimental study at 420C, 500 bars. M. S. thesis, Department of Earth Sciences, University of Minnesota.

Schmidt MW, Poli S. 1998. Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation. *Earth Planet Sci Lett* 163: 361-37.

Schwartz S, Allemand P, Guillot S. 2001. Numerical model of the effect of serpentinites on the exhumation of eclogitic rocks: Insights from the Monviso ophiolitic massif (Western Alps). *Tectonophysics* 342: 193-06.

Schwartz, S., Guillot, S., Reynard, B., Lafay, R., Debret, B., Nicollet, C., Lanari, P., & Auzende, A. L. (2013). Pressure-temperature estimates of the lizardite/antigorite transition in high pressure serpentinites. *Lithos*, 178: 197-210.

Sdrolias M, Müller RD. 2006. Controls on back-arc basin formation. *Geochem Geophys Geosyst* 7.

Turcotte DL, Schubert G. 2002. *Geodynamics*. Cambridge University Press.

Van Keken PE, Hacker BR, Syracuse EM, Abers GA. 2011. Subduction factory: 4. Depth-dependent flux of H₂ O from subducting slabs worldwide. *J Geophys Res* 116: B01401.

van Hunen J, van den Berg AP, Vlaar NJ. 2000. A thermo-mechanical model of horizontal subduction below an overriding plate. *Earth Planet Sci Lett* 182: 157-69.

van Hunen J, van den Berg AP, Vlaar NJ. 2004. Various mechanisms to induce present-day shallow flat subduction and implications for the younger Earth: a numerical parameter study. *Phys Earth Planet Inter* 146: 179-94.

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

Wagner LS, Anderson ML, Jackson JM, Beck SL, Zandt G. 2008. Seismic evidence for orthopyroxene enrichment in the continental lithosphere. *Geology* 36: 935.

Wang, H., Huismans, R. S., & Rondenay, S. (2019). Water Migration in the Subduction Mantle Wedge: A Two-Phase Flow Approach. *Journal of Geophysical Research: Solid Earth*, 124: 9208-9225.
https://doi.org/10.1029/2018JB017097

Watremez L, Burov E, d'Acremont E, Leroy S, Huet B, Le Pourhiet L, Bellahsen N. 2013. Buoyancy and localizing properties of continental mantle lithosphere: Insights from thermomechanical models of the eastern Gulf of Aden. *Geochem Geophys Geosyst* 14: 2800-817.

Wilson, C. R., Spiegelman, M., Van Keken, P. E., & Hacker, B. R. (2014). Fluid flow in subduction zones : The role of solid rheology and compaction pressure. *Earth and Planetary Science Letters*, 401: 261-274.
https://doi.org/10.1016/j.epsl.2014.05.052.

Wu S, Wang Y, Qian X, Asis JB, Lu X, Zhang Y, Gan C. 2022. Discovery of the Late Cretaceous Barru adakite in SW Sulawesi and slab break-off beneath the Central Indonesian Accretionary Complex. *J Asian Earth Sci* 232: 105214.

Zwaan F, Brune S, Glerum A, Vasey DA, Naliboff JB, Manatschal G, Gaucher EC. 2023. *Rift-inversion orogens are potential hotspots for natural H2 generation*. https://doi.org/10.21203/rs.3.rs-3367317/v1

**Cite this article as:** Gauthier A, Larvet T, Pourhiet LL, Moretti I. 2024. Water budget in flat vs. steep subduction: implication for volcanism and potential for H₂ production, *BSGF - Earth Sciences Bulletin* 195: 26.