Mantle serpentinization and associated hydrogen flux at North Atlantic magma-poor rifted margins

Hydration of ferrous-iron–bearing minerals in mantle peridotites, a process known as serpentinization, is one of the main mechanisms by which H₂ is produced in nature (e.g., Klein et al., 2020). Serpentinization has drawn much attention regarding its ability to provide the metabolic energy for subsurface chemosynthetic microbial communities (e.g., Ménez et al., 2012), thereby driving the emergence of early Earth ecosystems before the advent of photosynthesis (e.g., Canfield et al., 2006), and its potential to act as a sustainable, carbon-free energy resource (e.g., Lefeuvre et al., 2021, 2022).

Wide-angle seismic data and International Ocean Discovery Program (IODP) drilling at magma-poor rifted margins of the North Atlantic Ocean suggest that the mantle is exhumed at 100–180-km-wide continent-ocean transition zones (COTs; Fig. 1). At the West Iberia (WIM) and Newfoundland (NF) conjugate margins, IODP samples consist of hydrated peridotite with only minor magmatic intrusions (Jagoutz et al., 2007; Tucholke et al., 2007). Analysis of serpentinized peridotites indicates that serpentinite hosted hydrothermal systems that were active during extension (e.g., Klein et al., 2015). Mantle serpentinization is well studied at the rock sample scale and is thought to sustain chemosynthetic communities (Klein et al., 2015). However, the production of H₂ at the lithospheric scale is poorly constrained, and its global significance for H₂ flux during the breakup of Pangea has not been addressed (Merdith et al., 2020).

Laboratory experiments have shown that the serpentinization reaction is largely controlled by the reaction temperature and the water-rock ratio (McCollom et al., 2016). Brittle damage zones are the primary conduits for seawater circulation and therefore control the water-rock ratio (e.g., Caine et al., 1996; Pérez-Gussinyé and Reston, 2001). Thus, understanding the temperature and deformation history during the formation of COTs is key to estimating the H₂ flux. We present a two-dimensional (2-D) thermomechanical model that simulates the tectonic structure of the WIM-NF margins, and we coupled this model with the chemical processes of olivine serpentinization to evaluate how deformation, serpentinization, and hydrogen production occur during mantle exhumation. We estimated the H₂ flux rate at the COTs of this margin pair and showed its dependency on extension velocity and hydrothermal cooling. We then calculated the area of the COTs along the North Atlantic magma-poor rifted margins from seismic profiles (Fig. 1) and estimated the possible accumulated H₂ production during the opening of this sector of the North Atlantic Ocean.

To approximate the temperature and deformation during rifting, we used a thermomechanical model (Pérez-Gussinyé et al., 2020), further described in the Supplemental Material¹. Mantle serpentinization was modeled by assuming that the rate of serpentinized rock production, , is a function of temperature, T, as

where ε̇_p is plastic strain rate, ε_p is plastic strain, and ε̇_p^cri and ε_p^cri are the strain rate and strain criteria for serpentinization. Following Emmanuel and Berkowitz (2006), A was set to be 10^–10 s^–1, b_s was set to be 2.5 × 10^–4 °C^–2, and c_s was set to be 270 °C. We also tested another kinetic function given by Malvoisin et al. (2012). The difference between the two was small and did not affect our main results (Fig. S2). We assumed that the supply of water by fluid flow is more rapid than the uptake of water by serpentinization, as there is no evidence for serpentinization under water-limited conditions (Barnes and Sharp, 2006).

By setting the plastic strain rate and plastic strain criteria for serpentinization, we assumed that only active faults act as effective conduits for fluid circulation. Plastic strain is given by a non-associated plastic flow law (see the Supplemental Material), which allows for the examination of the yield surface according to the Drucker-Prager yield equation (Equation S6) and for deformation localization in shear zones or “faults.” The value ε̇_p^cri was set to be 10^-14 s^–1, while ε_p^cri values of 0.1 and 0.5 were explored. The reason for this choice lies in the fact that fluids tend to flow toward active faults due to negative tectonic pressure (e.g., Mancktelow, 2008), but it is controversial how fluid flow is distributed in deformation zones due to the complexity of fault permeability (e.g., Caine et al., 1996). The values of 0.1 and 0.5 were set to approximate the cases of active faults acting either as distributed or localized conduits for fluid flow, respectively.

For a given numerical time step, the rate of mantle serpentinization along the modeled section, M_t, is

where is the serpentinization rate for a given element i, S_i is the area of an element, ρ_m is the density of the mantle, and the summation is done over all elements, nel. As our model is 2-D and runs in the direction of extension, M_t is in units of kg/(yr × km), where km refers to the direction perpendicular to extension.

We first show how the mantle is deformed and exhumed in space and time in a reference case; we then show the plausible H₂ flux during mantle exhumation. It is important to note that the model we present, although it fits the available observations well, may not be unique. The model parameters that affect crustal thinning, such as the crustal rheology and initial Moho temperature, are discussed and justified in the Supplemental Material.

Figure 2C shows the model evolution for an initial Moho temperature of 600 °C and a full extension rate of 15 mm/yr. After 30 km of extension, the deformation localizes into four main faults (C1.1 in Fig. 2C). As extension increases, one of the faults becomes dominant. During slip of the active fault, the footwall gradually cools and increases in strength, while the hanging wall gradually heats up and decreases in strength. This process causes the deformation to migrate oceanward in a series of younger sequential faults. The occurrence of this deformation pattern and its duration depend on parameters such as the initial rheological profile and temperature gradient (Brune et al., 2014; Liu et al., 2022). The deformation pattern explains the tectonic asymmetry observed at the WIM-NF conjugate margins and the younging of the synrift sequences at the WIM (Lymer et al., 2019; Ranero and Pérez-Gussinyé, 2010).

Early serpentinization occurs in the footwall of active faults (C2.1 in Fig. 2C). This is because the mantle material in the fault’s footwall is carried to shallow depths during exhumation, where the temperature is cool enough for serpentinization. Given the low friction strength of the serpentinite, its occurrence localizes plastic deformation and forms a long-lived detachment fault (C2.2 in Fig. 2C). As extension increases, the detachment footwall back-rotates and is offset by a series of small landward-dipping conjugate faults (C2.3 in Fig. 2C). When the detachment becomes inactive, the model starts to exhume mantle symmetrically, with faults dipping in alternate directions (C3.1 in Fig. 2C). Those crosscutting faults are similar to the flip-flop detachment faults previously described (Reston and McDermott, 2011; Theunissen and Huismans, 2022) where the newly formed fault cuts the old fault plane in two.

The average rate of mantle serpentinization, M_t, is 0.33 × 10⁹ kg/(yŕ km) (Fig. 2A). The M_t values fluctuate; they are high when a new fault emerges, and fresh mantle rocks are hydrated, and they are low in the late stage of fault activity, when mantle rocks satisfying the conditions for serpentinization have been completely hydrated.

Mapping the mass of mantle rock into the amount of produced H₂ was done by considering two end-member olivine serpentinization reactions:

1. all Fe partitions into magnetite:

   
2. all Fe ends up in the Fe-rich serpentine end-member hisingerite (Fe₂⁺³Si₂O₅(OH)₄) or cronstedtite ((Mg,Fe⁺²)₂Fe₂⁺³SiO₅(OH)₄):

   

   

Assuming the above chemical processes contribute equally to the production of H₂, 15 moles of olivine are consumed to produce 1 mol of hydrogen; i.e., H₂ production per kilogram of rock is 454 mmol. Those equations will give upper-limit values because any partitioning of Fe ^{+ 2} into serpentine or brucite will lower the H₂ yield. Previous studies have suggested that about half of the Fe in serpentine is trivalent and contributes to H₂ production (Klein et al., 2009; Andreani et al., 2013). With this assumption, the amount of H₂ that can be produced from 1 kg of mantle rock is 227 mmol. This value is consistent with the range of 150–300 mmol/kg obtained from thermodynamic modeling of serpentinized mantle rocks from the WIM (Albers et al., 2021). Combining this estimation with the M_t (Fig. 2A), we predict a H₂ production rate, H_2t, of (5–10) × 10⁷ mol/(yr × km). This value is of the same order of magnitude as previous models for ultraslow-spreading centers (Cannat et al., 2010).

Figure 3 shows the effect of extension velocity and hydrothermal cooling. As the extension rate increases, the speed of mantle exhumation increases, leading to an increase in the H_2t. The role of hydrothermal cooling is more subtle. H_2t increases with increasing hydrothermal cooling in general. However, the effect is more pronounced for the highest extension velocity case, because the effect of cooling on the isotherms is more noticeable (Fig. 3). It is also interesting to note that the magnitude of the oscillation in H_2t, denoted by error bars, is smaller in models without hydrothermal cooling effects (C.1 and C.2 in Fig. 3C), due to the relatively thinner brittle lithosphere in these models and the smaller offset and shorter life span of faults (Fig. 3B).

Our model explains a range of observations at the WIM-NF margins (Fig. 4):

1. The model forms a subhorizontal deformation remnant at the distal margin, which is overlain by several faulted blocks. This structure is like the S reflectors observed on the seismic profile (e.g., Dean et al., 2015).

2. The spacing and uplift magnitude of the peridotite ridges in the model are consistent with observations (Fig. 4).

3. The extent of serpentinization within the basement in our reference model is 6–8 km, which agrees well with the most recent wideangle models published along the WIM (e.g., Grevemeyer et al., 2022) and also those in ultraslow-spreading environments, such as the Southwest Indian Ridge (Cannat et al., 2019).

Given that our model reasonably fits the available observations, we believe that it can provide a plausible paleotemperature and deformation history and therefore potential serpentinization and associated H_2t.

The results allowed us to estimate cumulative H₂ emission during the transition from rifting to spreading in the southern North Atlantic Ocean. The reference model takes 12 m.y. to reproduce a structure like the COT at the WIM-NF conjugate margin, which has a width of L = 180 km (Fig. 4). During this time, the accumulated hydrogen production, H_{2_acc} = H_2t× 12 m.y., is 9.75 × 10 ¹⁴ mol per kilometer of rift. To generalize our results to other North Atlantic magma-poor continental margins with similar rifting backgrounds, we needed to estimate the area of the COTs in the North Atlantic. We integrated published wide-angle seismic profiles to determine the width of the COTs in some regions and obtained constraints between profiles using interpolation (Fig. 1A; see the Supplemental Material). In this way, we obtained the area of COTs in the North Atlantic of S = 8 × 10⁵ km². Thus, the estimated accumulated H₂ flux from mantle serpentinization prior to the opening of the North Atlantic Ocean is

A fraction of this H₂ would have been consumed by reduction of CO₂ to hydrocarbons (Proskurowski et al., 2008), mostly CH₄, for which 4 moles of H₂ are required per mol of CH₄ produced. Previous studies have indicated molar H₂ to methane emission ratios related to the hydration of mantle rocks between 10:1 and 1:1 (Cannat et al., 2010). A 10:1 ratio would imply an H₂ flux of 3.1 × 10¹⁸ mol and a CH₄ flux of 3.1 × 10¹⁷ mol. A 1:1 ratio would equate to a flux of both compounds of 8.6 × 10¹⁷ mol. Hence, the amount of methane released from mantle serpentinization during the opening of the North Atlantic Ocean may have approached 10¹⁸ mol. We consider it unlikely that sulfate reduction was a significant throttle on H₂ release, since it should have been present only in very small amounts at the temperature of serpentinization.

Despite the uncertainties introduced by the treatment of serpentinization and the related release of reducing power, our results indicate very clearly that the COTs are sites of significant H₂ and methane production. Considering that the formation of COTs in the North Atlantic took ~30 m.y. (Barnett-Moore et al., 2018), we obtained a flux of H₂ as high as 1.4 × 10¹⁷ mol/m.y. This would indicate that the H₂ release rate was roughly 25% of what it is today for the global system of oceanic spreading ridges (7 × 10¹⁷ moles/m.y.; Merdith et al., 2020). As to the fate of the H₂ released by serpentinization in COTs, we expect that most of it was utilized by microorganisms in the subseafloor (e.g., Klein et al., 2015; Worman et al., 2020). Hydrogenotrophic sulfate reduction and methanogenesis are the most likely anaerobic catabolic reactions. Any venting of hydrothermal solutions would have added H₂ and CH₄ to the water column, which would have boosted chemolithoautotrophic biomass production there. The considerable surplus in the geogenic reducing power during COT formation is likely to have prominently enhanced microbial activity in, at, and above the seafloor within the region of the opening North Atlantic, but also at a global scale. All of this suggests that the production of H₂ and CH₄ during COT formation is an important part of the global ocean biogeochemical cycle, but one that has not been considered in current global flux estimation models (Merdith et al., 2020). Additionally, at COTs, unlike at mid-ocean ridges, H₂ may have accumulated in sealed reservoirs within the sedimentary layers, if they were within the temperature range where H₂ is relatively inert (100–200 °C; Lefeuvre et al., 2022). In the future, more detailed simulations including hydrogen production, migration, and sedimentation history are needed to analyze the potential of COTs as a source of green energy.

This work was supported by German Research Foundation grant 396827560. Z.L benefited from the Program for JLU (Jilin University, China) Science and Technology Innovative Research Team (No. 2021TD-05). We thank the editor Gerald Dickens and anonymous reviewers for their insightful feedback.