Field-Based Evidence for Intra-Slab High-Permeability Channel Formation at Eclogite-Facies Conditions During Subduction.

Fluid release from subducting oceanic lithosphere is a key process for subduction zone geodynamics, from controlling arc volcanism to seismicity and tectonic exhumation. However, many fundamental details of fluid composition, flow pathways, and reactivity with slab-forming rocks remain to be thoroughly understood. In this study we investigate a multi-kilometer-long, high-pressure metasomatic system preserved in the lawsonite-eclogite metamorphic unit of Alpine Corsica, France. The fluid-mediated process was localized along a major intra-slab interface, which is the contact between basement and cover unit. Two distinct metasomatic stages are identified and discussed. We show that these two stages resulted from the infiltration of deep fluids that were derived from the same source and had the same slab-parallel, updip flow direction. By mass balance analysis, we quantify metasomatic mass changes along this fluid pathway and the time-integrated fluid fluxes responsible for them. In addition, we also assess carbon fluxes associated with these metasomatic events. The magnitude of the estimated fluid fluxes (104-105) indicates that major intra-slab interfaces such as lithological boundaries acted as fluid channels facilitating episodic pulses of fluid flow. We also show that when fluids are channelized, high time-integrated fluid fluxes lead to carbon fluxes several orders of magnitude higher than carbon fluxes generated by local dehydration reactions. Given the size and geologic features of the investigated metasomatic system, we propose that it represents the first reported natural analogue of the so-called high permeability channels predicted by numerical simulations.

Introduction

Fluids produced by metamorphic devolatilization reactions exert major controls on mass transfer in the crust and mantle, as well as on the generation of arc magmas and new continental crust (Fyfe et al., 1978; Taira et al., 1998; Tatsumi & Eggins, 1995). However, even though element solubility in metamorphic fluids has been demonstrated by the examination of metamorphic veins and their associated alteration selvages (Ague, 1994a, 1994b, 2017; Bebout & Barton, 1993; Beinlich et al., 2010; Breeding & Ague, 2002; Gao et al., 2007; Masters & Ague, 2005; Penniston‐Dorland & Ferry, 2008; Spandler et al., 2003; Spandler et al., 2008; van der Straaten et al., 2012), how fluids escape from the subducting slab has been a longstanding issue and the fluid pathways remain vaguely understood (Connolly & Podladchikov, 2004; Faccenda, 2014; Faccenda et al., 2012; Hacker et al., 2003; Plümper et al., 2017; Wilson et al., 2014). While fluid production is due to metamorphic reactions (e.g., Schmidt & Poli, 2014), fluid expulsion is a mechanical process controlled by both the rheological properties of the solid matrix and by the hydraulic properties regulating the flow of the fluid through the matrix (Connolly, 2010). Metamorphic devolatilization reactions involve changes in density, porosity, permeability and fluid pressure which influence fluid flow (Connolly, 2010; Malvoisin et al., 2015; Plümper et al., 2017; Tian & Ague, 2014; Tian et al., 2018). Porosity formation is critical to the dynamics of fluid expulsion, and high fluid fluxes are expected in high‐permeability regions (e.g., regions developed in compacting two‐phase flow systems and lithological contacts; Ague, 1994a, 1994b, 2014; Connolly & Podladchikov, 2007; Oliver, 1996) that may evolve into fluid conduits (Ague, 2017; Connolly, 2010). Numerical studies have shown that compaction pressure gradients can divert the direction of slab‐derived fluids. In particular, compaction has the effect of enhancing updip fluid migration and localizing fluid fluxes along so‐called high‐permeability layers or compaction channels (Faccenda et al., 2012; Wilson et al., 2014). Such fluid channelization is a viable method to allow fluids to travel long distance (km‐scale) and transfer water and carbon toward the mantle wedge (Tian et al., 2019; Wilson et al., 2014). Nevertheless, natural case studies depicting the chemical‐mechanical feedbacks of fluid flow along coherent high‐pressure (HP) metamorphic units (i.e., blueschist and eclogite‐facies) are still limited (Giuntoli et al., 2020; Locatelli et al., 2018; Maggi et al., 2014; Plümper et al., 2017; Taetz et al., 2018) and most of the documented metasomatic rocks belong to mélange terranes (Bebout, 2007; Bebout & Barton, 1993; Breeding et al., 2004; Johnson & Harlow, 1999; Miller et al., 2009; Penniston‐Dorland et al., 2014; Sorensen et al., 1997; Spandler et al., 2008).

In this study, we investigate a regional scale metasomatic system (Martin et al., 2011; Piccoli et al., 2018; Piccoli et al., 2016; Vitale Brovarone et al., 2014) localized along a continuous lithological boundary (>13 km along strike) within an extremely well‐preserved example of coherent subducted oceanic to transitional Tethyan lithosphere. The investigated lithological boundary represented the onlap unconformity surface at the base of the Mesozoic post‐rift sedimentary cover rocks (basal unconformity, hereafter) in the HP unit of Alpine Corsica (Beltrando et al., 2014; Vitale Brovarone, Groppo et al., 2011; Vitale Brovarone, Beltrando, et al., 2011). We integrate previously published field and petrological data (Martin et al., 2011, 2014; Piccoli et al., 2016; Vitale Brovarone et al., 2014; Vitale Brovarone & Beyssac, 2014; Vitale Brovarone et al., 2018 ) with new geochemical data and new mass balance analyses to understand the source, direction, and mechanism of fluid flow along this surface. We conclude that the regional scale metasomatic front preserved in Alpine Corsica represents a fossil high‐permeability channel where fluids flowed with a slab‐parallel, updip direction via solid matrix deformation.

Geological Setting and Previous Results on HP Metasomatism in Alpine Corsica

The metasomatic system in this study is located in the lawsonite eclogite‐facies Monte San Petrone unit. This tectono‐metamorphic unit consists of a coherent section of the subducted Mesozoic slow‐spreading oceanic and transitional lithosphere of the Tethys Ocean (Beltrando et al., 2014; Vitale Brovarone, Groppo et al., 2011; Vitale Brovarone, Beltrando, et al., 2011) (Figure 1a) that underwent peak metamorphic conditions of about 500°C–550°C and ∼2.2–2.4 GPa (Caron & Péquignot, 1986; Pequignot, 1984; Vitale Brovarone et al., 2011) in the late Eocene (∼34 Myr; U‐Pb zircon, Martin et al., 2011; Lu‐Hf garnet, Vitale Brovarone & Herwartz, 2013). The internal structure of the Monte San Petrone unit consists of a basal body of serpentinite (hydrated mantle rocks) overlain by a laterally variable lithostratigraphy comprising (Figure 1b): i) slivers of continental basement rocks; ii) pillow metabasalts; and iii) Mesozoic metasedimentary rocks. Based on structural, lithostratigraphic, and metamorphic features, the contact between the Mesozoic metasedimentary rocks and the structurally underlying serpentinites has been interpreted as the basal onlap unconformity of the subducted Tethyan lithosphere, as also documented along several sections of the Alpine belt (Figure 1b; Beltrando et al., 2014; Vitale Brovarone, Groppo et al., 2011; Vitale Brovarone, Beltrando, et al., 2011).

Figure 1

(a) Sketch of the geodynamic context showing the subduction of transitional oceanic lithosphere (i.e. equivalent for the Monte San Petrone unit). Compaction and shear deflect fluid released from serpentinites up‐slope. Fluids are then channelized along subducted basal unconformity surfaces. Stars indicate the location of metasomatic rocks, both for Stage#1 and Stage#2. Figure modified after Angiboust et al. (2012). (b) Geological map of the Monte San Petrone unit showing the investigated basal unconformity surface (red line, dashed red line where the contact is poorly exposed) (modified after Pequignot, 1984; Vitale Brovarone, Groppo et al., 2011; Vitale Brovarone, Beltrando, et al., 2011).

Previous works have shown that rocks overlying the basal unconformity recorded intense metasomatism during prograde to peak metamorphism (Martin et al., 2011, 2014; Piccoli, 2017; Piccoli et al., 2018, 2016; Vitale Brovarone et al., 2014; Vitale Brovarone & Beyssac, 2014). A rather continuous, regional‐scale metasomatic front can be followed over about ∼13 km (Figures 1 and 2) (Piccoli et al., 2018, 2016; Vitale Brovarone et al., 2014). This metasomatic front overprinted the cover unit (thin continental crust slivers and metasedimentary rocks, Figure 2), whereas the serpentinized basement is unaltered (Piccoli et al., 2018). Two prograde stages of metasomatism consist of: i) a first, prograde, metasomatic stage characterized by the replacement of the rocks structurally overlying the basal unconformity with calc‐silicate rocks composed of diopside‐ and/or lawsonite‐rich assemblages (Stage#1 metasomatic rocks, Figure 2b; Martin et al., 2011; Piccoli et al., 2018; Vitale Brovarone et al., 2014); ii) a second metasomatic stage structurally and mineralogically superimposed on the previous one and characterized by intense carbonation (Figure 2c; Stage#2 metasomatic rocks in Piccoli et al., 2018, 2016). Here we summarize the main geochemical features of Stage#1 and Stage#2 rocks that are relevant for this study.

Figure 2

(a) Geological cross section along the A‐A′ profile of Figure 1b and field image of unaltered Mesozoic metasedimentary rocks. (b) Stage#1 metasomatism along the basal unconformity, field image and photomicrograph of the Stage#1 metasomatic rock (c) Stage#2 metasomatism superimposing on Stage#1 metasomatic rocks, field image showing the carbonation front replacing Stage#1 rock, and photomicrograph showing a replacement texture of Stage#1 assemblage by carbonate (Stage#2).

Stage#1 rocks formed at the expense of both Mesozoic metasedimentary rocks and thin continental crust slivers (Figure 2b). The metasomatic product consists of Ca, Al rich rocks (lawsonite and diopside rich) with no relict carbonate (precursor carbonate minerals were completely dissolved). Metasomatic lawsonite is characterized by complex trace element zoning in rare earth elements (REE), large ion lithophile elements (LILE), Cr, and Ti (Vitale Brovarone et al., 2014). On the basis of lawsonite geochemistry and garnet oxygen isotope composition, Martin et al. (2014) and Vitale Brovarone et al. (2014) propose that serpentinites were the main fluid source responsible for metasomatism, with contributions from Ca‐rich lithologies such as mafic or metasedimentary rocks.

Stage#2 metasomatism (i.e., carbonated metasomatic rocks) is superimposed on Stage#1 assemblages and is characterized by the precipitation of abundant Ca‐carbonate (calcite and locally preserved fresh aragonite), omphacite, and garnet (±lawsonite) (Figure 2c). In these rocks, microstructural relationships give corrosion‐like textures, notably with heterogeneous consumption and segmentation of diopside‐rich domains (relict of Stage#1), overgrown by omphacite. Ca‐carbonate + omphacite, Ca‐carbonate + lawsonite, and Ca‐carbonate crack‐seals and hydraulic breccias sealed by carbonate are also observed. The characteristic oxygen isotopic composition of carbonates indicates equilibration with mafic/ultramafic rock‐derived fluid (light O isotope values; δ18O around 10‰–12‰), whereas carbon isotopes indicate that C was more likely derived from metasedimentary carbonates (δ13C around 0‰) (Piccoli et al., 2016). This finding, together with mineral chemistry, suggests that, in analogy with the Stage#1 metasomatic event, fluids responsible for Stage#2 metasomatic events also had a mixed ultramafic/metasedimentary source. In the proposed model, fluids produced by major dehydration reactions (e.g. antigorite‐out reaction) in the deeper part of the slab (10–30 km below) percolate parallel to the slab and chemically evolve by reactive flow, becoming enriched in constituents such as Ca and CO2 (Figure 1a and Piccoli et al., 2018).

Methods

Sampling Strategy

Samples for this study were selected from the suite studied by Vitale Brovarone et al. (2014) and Piccoli et al. (2018), on the basis of their macroscopic and petrographic properties. Samples described in this work were collected from localities where metasomatic rocks are juxtaposed with serpentinite and were categorized as (i) protolith rocks, with no evidence of fluid‐rock interactions (metasedimentary rocks), (ii) Stage#1 rocks (carbonate free, lawsonite + diopside metasomatic rock) and (iii) Stage#2 carbonated metasomatic rocks. In detail: four samples of Mesozoic metasedimentary rocks, two samples of Stage#1 metasomatic rocks, and two samples Stage#2 carbonated metasomatic rocks were analyzed for bulk rock composition and used for mass balance analysis (supplementary information and Table S1, Piccoli et al., 2020 DATASET).

Mass Balance and Time‐Integrated Fluid Fluxes

Geochemical Reference Frame

Mass balance analysis is done in terms of constant mass of an immobile element in order to account for volume changes (Ague, 1994a, 1994b; Ague & Van Haren, 1996). Here we consider zirconium (Zr) as the immobile element and we use the weighted reference element procedure and statistical methods for mass balance calculation as reported by Ague (2011). The choice of Zr as reference immobile species is based on the extremely low Zr solubility in water‐rich HP fluids (Ayers & Watson, 1991; Breeding et al., 2004). Solubility of high field strength elements (HFSE) such as Hf and U is typically very low in HP metamorphic fluids; as shown in the following, their geochemical behavior is identical, within error, to that of Zr, thus further supporting this choice (Figure 3 and Figure S2).

Figure 3

Mass change in percentage for Stage#1 (a) and Stage#2 (b) metasomatic rocks relative to their precursors. Values for mass gain higher than 300% are reported in the box. In both cases, Zr is used as geochemical reference species. Only mass changes equal or larger than 50% are to be considered statistically significant.

Because the second metasomatic stage is superimposed on the first one, we subdivided the mass balance analysis into two steps:

The first step is to compare mass and volume changes between metasedimentary rocks (precursor rock) and Stage#1 metasomatic rocks (altered rock)

Then we compare Stage#1 metasomatic rocks (precursor rock) with Stage#2 metasomatic rocks (altered rock)

In this work, we discuss only mass changes of elements that have natural variability in the precursor rocks below 50%. Only mass changes larger than 50% are considered statistically significant.

Basic Equations

The equations used here to quantify mass and volume changes follow Ague (1994a, 1994b), Ague and Van Haren (1996), Philpotts and Ague (2009), and references therein. The gain and loss of elements during metasomatism will result in changes in rock mass and volume. The overall or total rock mass change (T mass) calculated using a reference immobile species i is given by:

T mass is a fraction; the percentage mass change can be found by multiplying it by 100. The superscripts 0 and ' indicate the precursor and altered rocks respectively. If rock mass has been lost, the concentration of i will be greater in the altered rock than in the precursor rock (residual enrichment of i). Conversely, if rock mass has been gained, the concentration of i will be less in the altered rock than in the precursor rock (residual dilution of i).

The mass change of a mobile species j estimated using species i as a reference is: where c indicates the concentration of j and is the fractional mass change for j.

The overall rock volume strain (or volume change) using a reference species i is given: where ε is the volume strain and and are the precursor and metasomatic rock density, respectively. Here, we model rock densities at peak conditions (500°C and 2.2 GPa) using either the software Perple_X (version 6.7; Connolly, 2009, modeled assemblage is: carbonates, phengite, quartz and lawsonite, consistently with observations on natural samples, see Supplementary Information) with the internally consistent thermodynamic database of Holland and Powell (1998, revised 2002) or by estimating rock densities from mineral modal abundances (Table 1).

Table 1

Rock Density Estimation

Modeled bulk for metasediments		Minerals	Density [kg/m3]
SiO	38.4	Garnet	4075
Al2O3	3.0	Diopside	3354
Fe2O3	2.1	Lawsonite	3120
MgO	0.9	Aragonite	2926.4
CaO	54.8	Rocks
MnO	0.2	Metasedimentary rocks a	2779
K2O	0.5	Stage#1 rock b	3237
Tot.	100	Stage#2 rock b	3143

Note. Modeled metasedimentary rocks composition (pre‐metasomatism protolith) is the average composition of analyzed samples.

Rock density estimated by phase diagram section calculations for the reported bulk composition.

Rock density estimated on the base of mineral modal proportion. Stage#1 rock: 50% clinopyroxene and 50% lawsonite; Stage#2 rock: 55% aragonite, 35% clinopyroxene, 5% garnet and 5% lawosonite (mineral abundances from Piccoli et al., 2018).

Time‐Integrated Fluid Fluxes

The time‐integrated fluid flux is the total amount of fluid that passes across an area during a given time interval (Baumgartner & Ferry, 1991; Bickle, 1992; Ferry & Dipple, 1991; Skelton et al., 1995). When the flow is advection dominated (distance of the geochemical front >100 m; Ague, 2014), time‐integrated fluid fluxes (q ) can be approximated by (Ague, 2014): where L GF is the distance of geochemical front propagation (measured parallel to the flow direction), s the concentration of the element i in the input fluid, is the concentration of the element i after the reaction, and n is the number of moles of i produced or consumed per unit volume rock (mol m−3). Here we set to 0, that is to say that all i contained in the fluid has been passed to the rock during reaction. This approach gives the most conservative results. We use the C concentration in aqueous fluids () in equilibrium with mafic rocks given by Molina and Poli (2000): CO2 = 3.4 wt.% (i.e., X CO2 of 0.05 or 930.68 mol m−3). It is worth stressing that choosing a higher C concentration in the input fluid (0.05 instead of, for example, X CO2 of 0.0007, 46.4 mol m−3, from Kelemen & Manning, 2015) yields an even lower estimate of fluid flux (Figure 5). The length scale of the geochemical front (L GF) can be constrained if the mechanism of fluid flow and flow direction are known. In this case study, cross layer diffusion (that would produce L GF of the order of meter‐scale) can be excluded because there is no symmetrical alteration of the underlying ultramafic rocks. Cross layer advection (i.e., fluids advecting from the underlying serpentinite, which would also give L GF of the order of meter‐scale) can be also excluded because there is no field evidence of dehydration in the serpentinized basement.

Figure 5

(a) Time‐integrated fluid fluxes for different geodynamic settings (modified after Zack & John, 2007), time‐integrated water fluxes in subduction zones (water fluxes* as re‐calculated from data reported in Cannat et al., 2010; van Keken et al., 2011) and time‐integrated fluid fluxes from this study (intra‐slab permeability channel). The width of the bar for water fluxes and intra‐slab fluid fluxes corresponds to minimum and maximum timescales of interaction (102–105 yr). (b) Time‐integrated Carbon fluxes (expressed in mols per km of slab) associated with dissolution/carbonation along high permeability channels (blue bar) calculated for different infiltrating fluid carbon content (X CO2 of 0.7 · 10−3 and 0.05 from data of Kelemen & Manning, 2015; Molina & Poli, 2000, respectively) and time‐integrated fluxes estimated by thermodynamic modeling for minimum and maximum timescales of interaction (Gorman et al., 2006).

Previous work has shown that fluid flow was advection dominated and the direction of the flow was parallel to the contact between serpentinite and overlying cover units (Piccoli et al., 2018). Herein we provide further evidence that the fluid flow was updip (km‐length scale) and that the regional occurrence of metasomatic rocks can be used to constrain the length scale of the geochemical front. This order of magnitude for the length scale (km‐length scale) is also consistent with the inferred depth of fluid source (Piccoli et al., 2018). In order to test the sensitivity of the model to variable L GF and parameters, a minimum (2 km) and maximum (10 km) values for L GF were explored (Table S2). The minimum distance of 2 km (across‐strike) is based on the (folded) exposed length of the basal unconformity along the profile in Figure 2a, whereas 10 km is based on the extent of the field area in map view (along‐strike). Here, only results for 930.68 mol m−3 (X CO2 = 0.05, Molina & Poli, 2000) are presented. Results for  = 46.4 mol m−3 (X CO2 = 0.0007, Kelemen & Manning, 2015), which give much higher flux estimates, are presented in Table S2.

Results

Precursor Metasedimentary Rocks

Precursor metasedimentary rocks are the metamorphosed Mesozoic sedimentary cover of the Tethys ocean basin (Figure 2). Like their equivalents in the Western Alps, they are characterized by the alternation of impure marbles and more terrigenous layers (Busigny et al., 2003; Garofalo, 2012; Vitale Brovarone, Groppo et al., 2011; Vitale Brovarone, Beltrando, et al., 2011). Mineral occurrence in carbonate rich metasedimentary rocks includes Ca‐carbonate, quartz, phengite, and accessory chlorite, opaque minerals and carbonaceous material. Lawsonite pseudomorphs replaced by white mica and chlorite are also observed (Piccoli et al., 2018). Metasedimentary rocks have a slightly negative REE pattern, with a negative Eu* anomaly of 0.2 (Eu* = Eu/((Sm + Gd)/2) (Figure S1a). Compositional homogeneity of precursor metasedimentary samples was tested using a concentration ratio diagram (Figure S2a). All measured samples have very comparable major and trace element concentrations (concentration ratios close to 1), with the only exception of Mn and Ni that display larger variability. Most important for this study, Ca, CO2, Zr and high field strength element (HFSE) ratios are close to one and satisfy the fundamental requirements of quantitative mass balance (Ague, 1994a, 1994b).

Average Percentage Mass Changes

Here we report the average chemical alteration that produced Stage#1 and Stage#2 metasomatic rocks for major elements, LILE, and HFSE. Details of REE mass change per individual samples are provided in the Supplementary Information (Figure S2c, S2d, and S2e). Importantly, values reported in this section are absolute mass changes expressed in percentage, where zero means no mass change, hence differences in measured concentration between precursor and altered rock (Table S1) are due to dilution or residual enrichment effects and not to an actual mass transfer.

Stage#1 Metasomatic Rocks

Stage#1 metasomatic rocks have high REE concentrations, and a chondrite normalized trend with negative slope similar to those of metasedimentary rocks (Figure S1a). The Eu anomaly ranges from slightly negative (Eu* = 0.4) to slightly positive (Eu* = 1.1). Major and HFSE concentrations in different samples of Stage#1 rocks are similar, whereas REE concentrations display a large compositional range (Figures S1b, and S2b). However, Ca, CO2, Zr and high field strength element (HFSE) ratios are close to one; hence for our purposes the selected samples meet the requirements for quantitative mass balance (i.e., a geochemical reference frame that can be assessed when using these samples to evaluate Stage#2 metasomatic mass changes, see Section 4.2.2).

Results of mass balance analysis between precursor unaltered metasedimentary rocks and Stage#1 metasomatic rocks indicate that mass loss dominates mass gain, such that the total mass loss relative to precursors is about −70 % and the total volume loss is −75 % (Table 2). Lawsonite + diopside rocks formed from precursor Mesozoic metasedimentary rocks display a significant mass loss in Ca and CO2 (−85% and −95%, respectively, Figure 3a) as well as significant losses of Si, Mn, and LILE, such as K, Rb, Cs, and Sr (Figure 3a). HFSE did not undergo significant mass change (with the exception of Th and Nb, possibly enriched because of the crystallization of accessory phases such as apatite and/or titanite) and display similar geochemical behavior compared to Zr (the reference species).

Table 2

Evaluation of Total Mass Change (T) and Total Volume Strain (εi)

Samples	ρ (kg m−3)	[Zr ppm]	Ti	εi
Unaltered metasedimentary rocks	2780 a	32.45	–	–
Stage#1 metasomatic rocks	3220 b	116.9 c	−70%	−75%
Stage#2 metasomatic rocks	3140 b	3.83 c	210%	220%
Unaltered metasedimentary rocks (Syros)	2780 a	13.9	–	–
Metasomatic Ep + Gln rock (Syros)	3190 a	40.5	−66%	−70%

Lws, lawsonite; Di, diopside; Ep, epidote; Gln, glaucophane(εi).

Density extrapolated from phase diagram sections.

Density calculated by mineral modal abundances.

Average value.

Stage#2 Metasomatic Rocks

Major and trace element concentrations in selected Stage#2 samples fall in a narrow range. All samples display a REE pattern that has a negative Eu anomaly (Eu* = 0.2; Figure S1a). Absolute concentrations of REE are higher than in metasedimentary rocks, but lower than in Stage#1 rocks (Figures S1a and S1b).

For Stage#2 metasomatism mass gain dominates mass loss, such that the total mass gain relative to precursors is about 210% and the total volume gain is ∼220% (Table 2). On average, Stage#2 rocks display a significant Ca and LOI (loss on ignition in Figure 3) mass gain of 560% and 1270%, respectively. The minimum mass changes considering variability in the precursor rocks are still very large: 230% and 439% for Ca and LOI, respectively. The main volatile‐bearing phase in Stage#2 rocks is Ca‐carbonate. Thus, we consider LOI mass gain as a proxy for CO2 mass gain. In addition, when converted to moles, the CaO and CO2 mass gains is apparently stoichiometric (1:1 proportion), thus supporting the hypothesis that CO2 is the dominant LOI component. Major elements Si, Mg, and Al do not undergo significant mass change (Figure 3b). Lower concentrations of Ti correspond to a mass loss of −58% (Figure 3b). Large mass gain is observed for Fe, Na, and Mn (median value 130%, 620%, and 905%, respectively) and LILE such as Sr and Pb which are highly compatible in carbonates.

Discussion

Element Mobility During COH‐Fluid‐Rock Interactions

Element mobility during dehydration reactions and reactive fluid flow exert important controls on arc magma chemistry and the geochemical budget of the subducting lithosphere (Bebout, 1991, 2007; Bebout & Penniston‐Dorland, 2016; Breeding et al., 2004; Ishikawa & Nakamura, 1994; Marchesi et al., 2013; Peters et al., 2017; Spandler et al., 2003). Fluids liberated at blueschist‐eclogite conditions are expected to be aqueous, containing only a moderate amount of LILE, Sr, and Pb (Hermann et al., 2006), consistent with the observation of high LILE/REE ratios in volcanic arc rocks (e.g., Gill, 2012; McCulloch & Gamble, 1991).

Previous isotopic studies on metasomatic rocks from the Monte San Petrone unit (Martin et al., 2014; Piccoli et al., 2018, 2016) suggested that the fluid responsible for both metasomatic stages was equilibrated with ultramafic/mafic rocks. Further interaction by reactive fluid flow along the subducted basal unconformity was responsible for carbonate dissolution during Stage#1 (for carbonate rich lithologies such as calcschists) and Ca, C, and Na enrichment during Stage#2. Our mass balance analysis confirms this interpretation and adds important new information. The nearly 1:1 molar enrichment of Ca and C (Ca: 1.62 mols/100 g of rock; CO2: 1.65 mols/100g of rock) during the Stage#2 event points to carbonate dissolution as a major Ca‐C source (Ague & Nicolescu, 2014). This implies that an addition of externally derived Ca is required although Ca in the precursor rock (hosted in the diopside and lawsonite) could have been internally recycled to form aragonite, omphacite and garnet. Among major elements with metasedimentary affinity, Na, Fe, and Mn are strongly enriched, likely hosted by the newly formed omphacite and garnet in Stage#2 rocks (Piccoli et al., 2018). A metasediment‐derived component is also indicated by the strong enrichment in LILE, Sr, and Pb (as suggested for other settings, Breeding et al., 2004; Hermann et al., 2006), whereas U‐Th, REE, and HFSE are not fractionated (Figure 3b and Figure S2c). Interestingly, U and Th do not show significant mass changes, suggesting a geochemical decoupling of Pb from Th and U, and overall U retention in the subducting slab during both dehydration (i.e., serpentinite dehydration; Peters et al., 2017; Pettke et al., 2018; Piccoli et al., 2019) and HP metasomatism.

The overall mass changes of most elements show opposite trends of mass gain/loss between Stage#1 and Stage#2 (Figures 3 and 4a). However, it is also evident that for major elements such as Na, Fe, Al, and Mn, mass loss and gain do not balance each other (Figure 4a). This observation suggests that the Stage#2 event did not only derive from an “in situ” dissolution re‐precipitation mechanism, but required an additional external fluid source (as also indicated by Ca‐CO2 enrichment).

Figure 4

(a) Bar diagram showing mass change expressed as oxide weight percent for Stage#1 (relative to precursor metasedimentary rocks, blue bars) and Stage#2 (relative to precursor Stage#1 rocks, orange bars). The sum of the two episodes (“overall mass change”) is represented with yellow bars. (b) Volume strain for the two metasomatic stages represented as a function of time, where Stage#1 precedes Stage#2 (as deduced from field cross‐cutting relationships).

Mass balance analysis also reveals that the two metasomatic stages are associated with significant volume changes (Table 2). Interestingly, the overall volume change brings rocks that underwent both metasomatic stages almost back to the initial volume (Figure 4b).

Fluid Fluxes in Subduction Zones

The origins, flux magnitudes, and transport mechanisms of metamorphic fluids are still subject to large uncertainties, and estimated time‐integrated fluid fluxes for subduction zones range from near zero to more than 103 (Figure 5a and Ague, 2014; Zack & John, 2007). These variations may reflect the different methods used to estimate paleo fluid flow. For our estimation, we used a value of C concentration in fluids of 3.5 wt.% (i.e. X CO2 of 0.05, Molina & Poli, 2000). The q obtained for both the two stages (considering minimum and median CO2 mass change) are in the order of up to (Figure 5a and Table 3) which corresponds to a carbon flux spanning from mol C/m2 to mol C/m2, for Stage#1 and Stage#2 respectively. These estimated q values are 2–3 orders of magnitude larger than the estimated fluid fluxes at the top of the slab (102 , Zack & John, 2007, see Figure 5a). When fluids are channelized (e.g., in veins), time‐integrated fluid fluxes near the top of the slab at subarc depth can be higher, up to 104 (Ague, 2014; Zack & John, 2007). This estimate coincides with our minimum estimation and is one order of magnitude lower than the maximum q value obtained.

Table 3

Evaluation of CO2 fluxes Considering a C Concentration in Fluids of 3.5 wt% and an Arc Length of 1050 km

	Initial content (wt.%)	Mass change	C eq (mol m−3)	Cinf (mol m−3)	L (m)	qTI (m3 m−2)	C Flux (mol m−2)	CO2 flux (mol km−1)
Metasediments	32	–	–	–	–	–	–	–
Stage#1 rocks	6	−95%	0	931	10000	2E+05	8E+08	−8E+17
					2000	5E+04	2E+08	−2E+17
Stage#2 rocks	28	1270%	0	931	10000	6E+05	2E+09	+2E+18
	–	440%	–	–	2000	4E+04	1E+08	+2E+17
Metabasalts a								−8E+08
Thermodynamic modeling b								−1E+12

Negative sign for CO2 release, positive value for CO2 uptake.

Value from Vitale Brovarone et al. (2018).

Value from Gorman et al. (2006), considering 10^5 years of fluid‐rock interaction.

Supercritical fluids at fore‐arc to subarc depth are dominantly aqueous (total amount of dissolved solids <30 wt.%, Hermann et al., 2006). Hence, one possible approach to test the reliability of our results is by comparing the estimated time‐integrated fluid fluxes with estimated cumulative water fluxes in subduction zones (van Keken et al., 2011). In other words, is the amount of water released by dehydration reactions from forearc to subarc depth is large enough to justify the fluid fluxes obtained here by mass balance analysis? The caveat for this approach is that the rock record only provides insights on the cumulative final result, and mass balance analysis does not provide information about the timescale of the fluid‐rock interaction. Attempts to constrain the timescale of metasomatic processes have been based on isotope diffusion chronometry; estimates span from very short (a few hundred years, Beinlich et al., 2020; Camacho et al., 2005; John et al., 2012) to relatively longer timescales (103–105 years, Penniston‐Dorland et al., 2010; Van Haren et al., 1996; Young & Rumble, 1993). Here we will use a minimum (0.01 cm/yr; Penniston‐Dorland et al., 2010) and maximum (10 cm/yr; Beinlich et al., 2020) isotope diffusion front propagation rates (which is perpendicular to the flow direction) to bound the minimum and maximum duration of the observed metasomatism. Considering a 10 m thick metasomatic halo (perpendicular to the fluid flow direction, Piccoli et al., 2018), the timescale in our case ranges from 102 to 105 years.

It is now possible to compare our time‐integrated fluid fluxes with the estimated water fluxes in subduction zones. van Keken et al. (2011) estimated a water flux of 35.9 Tg/Myr/m for the slow and cold subduction of Calabria, which, in terms of subducting rate, age of the oceanic lithosphere and thermal regime can be used as modern analog for Alpine Tethys subduction. Taking the minimum and maximum timescales for fluid‐rock interaction, the subducting slab can release from 109 to 1012 m3 of water. In an analogous way, we used the data from Cannat et al. (2010) to estimate the volume of subducted serpentinite for a slow‐spreading ocean (serpentinites are the main water source in subducted slow‐spreading oceanic lithosphere, Ulmer & Trommsdorff, 1995; see supplementary information for details on calculation). The total water loss for serpentinite dehydration reaction ranges from 108 to 1011 m3 of water. In both cases, for pervasive fluid flow, the released volume of water would result in low fluid fluxes (e.g. 0 ≤ q TI ≤ 102 if a slab section of 1000 m height and 400 km length is considered, i.e. length of the Calabria arc, see Table 2 in van Keken et al., 2011). On the other hand, if fluids were channelized in km‐scale porosity channels (e.g. cylindrical geometry with 2 km diameter, Tian et al., 2018), the resultant fluid fluxes range between 102 and 105 . Our estimated time‐integrated fluid fluxes fall within the upper range of values for the second scenario. It is noteworthy that high fluid fluxes are predicted for fluid flow along porosity channels in the mantle wedge, and can be plausibly realized in nature if deep‐sourced fluid plumes are deflected updip (Wilson et al., 2014). However, modeled peak subarc fluxes are in the order of 10−3 (Wilson et al., 2014). Even assuming a very generous porosity of 10%, for the minimum timescale of 102 year, the estimated time‐integrated fluid fluxes (Table 2) return unrealistic pore velocities (average flow velocities) of which are higher than all known subsurface flow systems. Smaller porosities yield even higher values; for example, the range is for porosity = 0.1%. Therefore, the timescale of 102 yr is considered unrealistic for this case study, and 105 yr should be considered as a more likely timescale for fluid rock interactions. In addition, knowing that the eclogitic peak is dated at 34 Myr and the blueschist at 37 Myr (Vitale Brovarone & Herwartz, 2013), and that the Nappes Supérieures were emplaced between 37 and 40 My (biostratigraphic constrains), a maximum timescale for HP fluid rock interaction can be assessed to be of the order of yr.

Altogether, considering the regional scale of the metasomatic front, the identification of regional‐scale basal unconformities as loci for fluid channelization, and the elevated time‐integrated fluid fluxes that are directly comparable with subduction zone water fluxes, we propose that the investigated field area represents the first reported field observation of an intra‐slab high permeability channel as hypothesized by numerical models.

Implications for Carbon Fluxes

Sequestration of volatiles during reactive fluid flow may strongly impact volatile fluxes (Vitale Brovarone & Beyssac, 2014; Wada et al., 2012). Most previous studies have focused on carbon released in response to decarbonation/carbonate dissolution processes (Ague & Nicolescu, 2014; Gorman et al., 2006; Kerrick & Connolly, 2001; Vitale Brovarone et al., 2018). However, a debate exists on the magnitude of carbon fluxes. For example, studies on metasedimentary and metamafic rocks concluded that most of the subducted carbon is retained in the slab (Collins et al., 2015; Cook‐Kollars et al., 2014). Others claimed that almost all carbon is returned to the surface (Kelemen & Manning, 2015). Here we attempt to reconstruct carbon budgets for an evolving reactive fluid flow considering both dissolution/decarbonation and carbonation reactions (Stage#1 and #2) along permeability channels and compare it with fluxes generated by local fluid release.

The Monte San Petrone unit recorded peak P conditions of 2.2–2.3 GPa, corresponding to forearc depths (70–80 km). For cold modern subduction zones (e.g., Honshu, Japan), Gorman et al. (2006) estimated a peak CO2 flux of 11.9 · 106 mol · km−1 yr−1 at ∼70 km depth (forearc), a depth that is comparable to our case study. Considering a time scale for fluid‐rock interaction of 102–106 yr, it results in a flux of 109–1013 mols · km−1 (Figure 5b). To compare our results with those of Gorman et al. (2006), the CO2 fluxes estimated in this study are multiplied by 1050 km (Honshu arc length, Table 3). The resulting flux is 2 · 1017 up to 8 · 1017 mols of CO2 km−1 for Stage#1 (carbonate dissolution Fig. 5b), and 2 · 1017 up to 2 · 1018 mols of CO2 km‐1 for Stage#2 (carbonation, Figure 5b). In order to account for other possible carbon sources at the regional scale, our results are compared to the study of carbon release from metabasalts from Alpine Corsica. However, the estimated fluxes from metabasalt decarbonation are significantly lower (8 · 108 mols of CO2 km−1, Vitale Brovarone et al., 2018). Figure 5b shows that the flux needed to account for the observed CO2 uptake during Stage#2 is several orders of magnitude higher than fluxes estimated by thermodynamic modeling for devolatilization reactions (Gorman et al., 2006), so we speculate that high carbon fluxes are mainly focused along high permeability regions due to the nonwetting property of CO2.

However, it is important to recall that in the model of Gorman et al. (2006), the subducting slab is a section of fast‐spreading oceanic lithosphere, whereas Alpine Corsica represents a subducted hyper‐extended continental margin. As a consequence, the initial volatile content of both the metasedimentary pile and altered oceanic crust is different. Moreover, the location on a P‐T diagram of the major dehydration reactions is different. For example, previous study and our field observations attest to the lack of lawsonite breakdown at peak metamorphic conditions (Pequignot, 1984; Vitale Brovarone, Groppo et al., 2011; Vitale Brovarone, Beltrando, et al., 2011). This means that, for the Alpine Corsica case, the water‐rich fluid source must come from deeper portions of the slab, plausibly from the antigorite out reaction, as previously suggested (Piccoli et al., 2018).

Fluid Pulses Along an Intra‐slab High Permeability Channel: A P‐T‐t Journey

Focused fluid flow is necessary to explain observed intra‐slab metasomatism as well as flux‐melting in the mantle wedge and arc magmas. Wilson et al. (2014) demonstrate that when compaction pressure gradients are included in numerical models, fluids released in the deeper part of the slab (e.g. dehydrating serpentinized mantle) are deflected updip and production of non‐linear porosity waves modulates the flow of fluids in the mantle wedge. The study of Malvoisin et al. (2015) shows that similar fluid transport mechanism can possibly be realized also in intra‐slab settings. The switch in type of metasomatism can be explained by three different processes: i) the different metasomatic stages were associated with different sources of the infiltrating fluid; ii) the fluid source was always the same, but there was a change in the direction of the flow (up‐T and carbonate dissolution, down‐T and carbonate precipitation); iii) fluid source and direction of the flow were the same, but the metasomatic products changed with changing P‐T conditions during continued subduction (P‐T‐t evolution). The first hypothesis (i) can be ruled out by previous studies on oxygen isotope composition of metasomatic silicates and carbonates that demonstrate that the water source was, in both cases, from ultramafic rocks (Martin et al., 2011, 2014; Piccoli et al., 2016).

It has been shown by numerical simulations and field evidence that up‐T fluid flow can lead to decarbonation and carbonate dissolution reactions (Gorman et al., 2006; Vitale Brovarone et al., 2018). Up‐T fluid flow can also lead to both decarbonation and carbonation at lithological interfaces, but these two metasomatic processes should be separate in space and not superimposed (i.e., CO2 is removed from one reservoir and redistributed in the one on top, Tian et al., 2019). A change in the direction of the flow (hypothesis ii) from up‐T (upward toward the hotter mantle wedge) to down‐T (updip) could potentially explain the co‐existence, at the same structural level, of carbonate dissolution and carbonation. However, the oceanic lithosphere underlying the basal unconformity is composed of fully hydrated serpentinites, which have very low permeability as shown by the absence of carbonation (Piccoli et al., 2018) and by recent experimental studies (Ganzhorn et al., 2019). Therefore, we also exclude the second hypothesis and we propose that in both cases the fluid flow was down‐T, slab parallel and updip.

In the third hypothesis, the fluid source was the same, but fluid infiltration took place episodically, by pulses, at different P‐T‐t conditions. In our interpretation, this model best explains the natural observations. We propose that the arrival of a first fluid pulse (t1 in Figure 6) is responsible for the carbonate dissolution and lawsonite + diopside precipitation. Then, during continued subduction, a second pulse arrived (t1′). This time the infiltrated rock is a carbonate‐free, metasomatic rock of broadly “mafic” composition, a condition favorable for carbonation reactions to take place. The presence of hydraulic breccias and crack‐seals suggests that fluid infiltration was episodic (Piccoli et al., 2018). Hence, we cannot exclude that multiple episodes of Stage#1 and Stage#2 (i.e., arrival of multiple fluid pulses) happened during subduction. Importantly, this also implies that this major intra‐slab high permeability channel remained active during a potentially long period of time (hundreds of thousands to millions of years).

Figure 6

Porosity evolution during fluid flow along an intra‐slab permeability channel. At t1, the arrival of the reaction front of a first fluid pulse generates porosity. At t2, in analogy with the model of porous media flow of McKenzie (1984), the reacted volume compacts against the impermeable layer (in this case study, the serpentinite). Compaction causes overpressure that overcomes the tensile strength. In t1′ new porosity is formed by decompaction and vein network formation, which allows the arrival of a second fluid pulse.

Possible Mechanism for Fluid Flow

The physical mechanisms responsible for the formation of intra‐slab permeability channels have yet to be understood. The two stages of metasomatism in the Monte San Petrone unit are characterized by large mass and volume changes. We quantify volume changes during metasomatism by mass balance analysis (Equation 3). In Table 2 absolute mass and volume changes for the Stage#1 and Stage#2 metasomatic events are reported. For comparison, total mass and volume changes associated with carbonate dissolution from the subduction complex in Syros island (Greece) are also reported (Ague & Nicolescu, 2014). The results for carbonate dissolution give quite consistent results between Alpine Corsica and Syros, with a total volume changes of ca. −70%. Volume changes associated with rock carbonation in Alpine Corsica are extremely large, around +200%. However, increasing the volume by +200% after loss of −70%, returns the reacted volume to near the initial volume (Figure 4b). These data suggest that matrix deformation (t2 in Figure 6) might have favored fluid flow and played an important role in the formation of the regional scale metasomatism. Actually, the compaction after the passage of the fluid pulse(s) (t2 in Figure 6) creates a flow barrier that can divert the new upcoming fluid pulses. Since the rock record shows the arrival of a second fluid pulse or pulses, a mechanism capable of creating new porosity is required. Our field observations (e.g. crack‐seals, hydraulic breccias, Piccoli et al., 2018) lead us to speculate that viscous decompaction during Stage#1 metasomatism generates overpressures that exceed the tensile strength of rocks to damage them (t1, Figure 6, Connolly & Podladchikov, 2007). These damaged weak zones can be preserved (e.g., as unhealed grain boundaries) even during the subsequent compaction period (t2, Figure 6). The arrival (t1′, Figure 6) of the second fluid pulse can exploit this pre‐existing weak and deformable zone as a flow channel to cause the overprint of Stage#2 metasomatism. We thus propose that fluid infiltration associated with matrix deformation (compaction and decompaction) can be an efficient mechanism to channelize fluids, maintain high fluid fluxes, and sustain the reaction progress.

Conclusions

Metasomatic rocks preserved in HP‐LT metamorphic units are the ultimate evidence for the passage of large amounts of fluids through rocks during subduction. Updip fluid flow has important consequences for the generation of arc magmatism. Actually, numerical models show that updip fluid migration can generate higher volatile fluxes (H2O and CO2) at sub‐arc depths (Tian et al., 2019; Wilson et al., 2014). Mass balance analysis of metasomatic rocks formed along a slab‐scale lithological interface indicates that extremely high fluid fluxes were realized within the slab. In our interpretation, the investigated metasomatic system represents a first report of a natural analogue of an intra‐slab high permeability channel, possibly active over hundreds of thousands of years, consistent with numerical simulations. Moreover, we show that when fluids are channelized, high time‐integrated fluid fluxes lead to carbon fluxes several orders of magnitude higher than carbon fluxes generated by local dehydration reactions, hence favoring the return of carbon back to the surface. Nevertheless, reactive fluid flow along an evolving metasomatic system can lead to CO2‐rock interactions and carbonation, which can sequester as much carbon as the amount liberated from decarbonation/dissolution reactions. Both dissolution‐carbonation reactions have to be taken into account when attempting to assess the net carbon flux in the context of global carbon cycle.

Lastly, we speculate that matrix deformation (compaction and decompaction) could have been the mechanical process allowing fluid flow. It is tempting to compare the observed peak metamorphic compaction‐decompaction to the porosity wave fluid flow model. While flow patterns matching with the porosity wave model have been recognized in asthenospheric setting (Bouilhol et al., 2011), no field based study reports evidence for such a process in intra‐slab settings. However, the possibility exists that porosity waves can be generated at depth by antigorite breakdown reactions and propagate upward (Connolly, 1997; Malvoisin et al., 2015; Skarbek & Rempel, 2016). More studies are needed to better investigate this possibility. We conclude that models including both mechanical and geochemical aspects of fluid flow are needed to better understand geochemical cycling in subduction zones.

Supporting Information S1

Click here for additional data file.