Water budget in flat vs. steep subduction: implication for volcanism and potential for H₂ production

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.

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 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 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.

The pTatin2D code is a thermomechanical modelling software specifically designed for studying long-term geodynamic 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.

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).

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).

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 ⁠. 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

which depends on cohesion (C) and friction angle (ϕ) or the maximum plastic strength σ_max fixed at 400 MPa (Watremez et al., 2013), the viscosity is adjusted accordingly to

The temperature is obtained solving the heat conservation

for temperature (T), using the thermal diffusivity (κ), the heat capacity (c_P), 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.

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).

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.

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.

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.

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

where n represents the number of markers in the element, S denotes the surface area of the element [m²], and ρ_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_H2O < 0), some markers inside might already be full (bw > wk). In such 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_H2O > 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_H2O > 0).

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.

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 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 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:

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

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 decreases 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).

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 (Grevemeyer 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).

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 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).

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).

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_PA) (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.

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 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 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 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.

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 exclusively 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_PA around 200 km. They also highlighted the uncertainties associated with the percolation constant. While this constant is geologically poorly constrained, 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.

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.

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 sab, 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.

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 $VpVs$ 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.

The combustion of Dihydrogen (H₂) 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₂ and/or consume a lot of energy. In contrast, natural sources of H₂ (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₂ 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₂ (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₂ production is well known in the formation of low-grade serpentines (lizardite and chrysolite), H₂ 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₂ production are unlikely. However, in recent studies Boutier et al. (2021) and Vitale Brovarone et al. (2020), have observed fluid inclusions of H₂ 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₂ was produced in the antigorite stability field which implies deep production of H₂ and suggests that H₂ 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₂ 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₂. Therefore, it is likely that, just like in the case of serpentinization, these reactions, which are also linked to water reduction, generate H₂ through the oxidation of Fe²⁺ present in peridotites under HP-HT conditions.

In summary, recent studies seem to show that mantle wedge hydration reactions could generate H₂. There are currently no equations to calculate the amount of generated H₂ based on the amount of consumed water in the mantle wedge. However, we can assume that the more water is consumed, the more H₂ is produced and use the amount of consumed water as a proxy for H₂ 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₂ 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₂ production during low-grade serpentine formation which is a better thermodynamically constrained redox reaction.

Our results highlight the competition between H₂ 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₂ production in flat subduction. Additionally, in a flat subduction context, the quantity of H₂ 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₂ exploitation, as the limitation of hydration zone extent increases H₂ concentration. Nevertheless, note that this study focused on H₂ production at depth, but cannot predict surface releases, as we do not model the H₂ migration from the deep source to the surface.

Our simulations have highlighted the importance of sediments in mantle hydration and consequently in H₂ 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…).

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.