Introduction

Increase of pressure and temperature during prograde metamorphism along a subducting slab leads to a series of breakdown reactions involving hydrous phases (e.g. Schmidt and Poli 2003) and is associated with major variation in rock density, rock volume and rock porosity. This process liberates varying amounts of free fluids (dominantly aqueous fluids) responsible for element transfer between crust and mantle in subduction zones. Fluids in turns play an important role in arc magmatism, geochemical cycles, energy budgets, mantle composition and seismicity. Studying the scale of fluid migration and the degree of fluid-rock interaction at depth provides information to quantify the fluid budget in the subducting slab environment (Sorensen and Barton 1987; Bebout 1991a; Philippot and Selverstone 1991; Selverstone et al. 1992; Nadeau et al. 1993; Spandler et al. 2011).

It has been demonstrated that large amounts of aqueous fluids in subduction zone can be generated by dehydration of mafic oceanic crust and associated hydrated oceanic lithospheric mantle (Ito et al. 1983; Peacock 1990; Ulmer and Trommsdorff 1995; Poli and Schmidt 2002; Spandler et al. 2003; Angiboust and Agard 2010). Aqueous fluids released from the lower layers of the subducting slab necessarily infiltrate upper layers of the subducting oceanic plate (Zack and John 2007) and eventually infiltrate the overriding mantle wedge. Propagation of fluids at depth may occur by interconnected vein networks and channelization (Miller and Cartwright 2000; Hacker et al. 2003; Miller et al. 2003; Zack and John 2007; John et al. 2008; Angiboust et al. 2014; Taetz et al. 2016) or by pervasive fluid flow (Bebout 1991b; Bebout and Barton 1993; Konrad-Schmolke et al. 2011). Pervasive fluid flow is hereby defined as fluid migration around individual mineral grains through an interconnected porosity, opposed to channelized flow that implies preferential fluid motion in one of more high permeability conduits, such as fractures and lithological contacts (Ague 2014). The fluid flow mechanisms may have different dynamics in the downgoing slab, at the slab-mantle interface or in the hanging-wall mantle wedge (e.g. Konrad-Schmolke et al. 2011; Angiboust et al. 2020). A better understanding of the dynamics of fluid-rock interaction and the associated permeability is crucial to quantify fluid-driven mass transfer processes in subduction zones (Ingebritsen and Manning 1999, 2003, 2010; Manning and Ingebritsen 1999; Konrad-Schmolke et al. 2011; Ganzhorn et al. 2019).

Stable isotopes are particularly suited for investigating the evolution of aqueous fluids in metamorphic terranes (e.g. Baumgartner and Valley 2001). They have been successfully applied to constrain the origin of subduction fluids and to describe fluid-rock and fluid-rock-melt interaction in subducting slabs (e.g. Barnes et al. 2014; Barnicoat and Cartwright 1997; Bebout and Barton 1989; Putlitz et al. 2000; Sharp et al. 1993). Distinct rock types are each characterized by a limited range of bulk δ18O values, and thus the δ18O composition of the derived fluids (water) at equilibrium can be estimated at a given temperature. In refractive metamorphic minerals that are robust to successive re-equilibration, such as garnet, mineral-fluid interaction will results in intra-crystalline δ18O variations when the reactive fluid is out of equilibrium with the bulk rock. Such variations in mineral oxygen isotope composition are a key tool to evaluate fluid-rock interaction with externally derived fluids and its relative timing during pressure–temperature (PT) evolutions (e.g. Errico et al. 2013; Gauthiez-Putallaz et al. 2016; Martin et al. 2014; Nadeau et al. 1993; Page et al. 2014, 2019; Putlitz et al. 2000; Rubatto and Angiboust 2015; Russell et al. 2013). Advances in thermodynamic modelling and underpinning experiments provide a detailed framework on fluid production and mineral δ18O change in single rock types (Kohn 1993; Zheng 1993; Manning 2004; Kessel et al. 2005). Recent studies have enabled major improvement in numerical modelling of petrological and isotopic systems by coupling thermodynamic and δ18O simulations (Vho et al. 2019, 2020a). With this approach, the effect of fluid interaction between multiple rock types of contrasting composition can be quantitatively evaluated along a given PT metamorphic path, allowing monitoring of variables of interest such as fluid-rock ratios and δ18O isotopic variations at the rock and mineral scale (Vho et al. 2020b).

Common dehydration reactions in metamorphic environments generate fluid and grain-scale porosity (Connolly 2010). Because the volume of fluid and minerals produced by devolatilizing reactions is greater than the volume of the mineral reactants, this pressurizes the porosity of the rigid host rock (Ferry 1994). This change in volume forms a micro-porous network, which permeability relies on the hydraulic connectivity among fluid-filled porosity (Miller et al. 2003; Zack and John 2007; Konrad-Schmolke et al. 2011). Even if fluids are poorly drained in deep environments (Connolly 2010), pervasive fluid flow will eventually set up in response to buoyancy forces and gradients in fluid pressure, where the interconnected porosity is sufficient to allow fluid motion (Oliver 1996; Ague 2003). The related fluid flux strongly depends on the rate of fluid production, on how deformation accommodates the rock volume variation and on the drainage efficiency (Connolly 2010). Fluid gradients and associated permeabilities, which vary in time and space, must be effective across the length of the fluid-flow system (Oliver 1996; Ingebritsen and Manning 2010) and remain sustained over sufficient time to accommodate the variations recorded by the exhumed metamorphic rocks (Manning and Ingebritsen 1999), such as a shift in stable isotopes or variation in major and trace elements bulk composition. Thus, on short time scales, a transient rock permeability may reach values significantly in excess than those expected (Ingebritsen and Manning 2010) and initiation of fluid flow will start as a consequence of fluid pressure dissipation (Manning and Ingebritsen 1999).

The European Alps are a world-class laboratory to observe relicts of fluid-rock interaction processes, with occurrence of a broad variety of metamorphic rocks exhumed from a paleo subduction zone. We investigated an association of metasediments and mafic felses (massive metamorphic rocks lacking clear schistosity or foliation) from the Theodul Glacier Unit, embedded within the Zermatt-Saas metaophiolite in the Western Alps. Both tectonic units reached eclogitic facies metamorphic conditions during Alpine convergence at ~ 45–50 Ma (Duchêne et al. 1997; Rubatto et al. 1998; Lapen et al. 2003; Bovay et al. 2021). Garnet chemical zoning and garnet δ18O variation are used to infer the record of fluid-rock interaction along the TGU metamorphic PT path. Thermodynamic and δ18O isotope modelling provides quantitative information on fluid source, fluid-rock ratio, time-integrated fluid flux and rock rheological variation. Modelling the variation in rock volume allows investigation of the origin of fluid transfer initiation and of rock permeabilities.

Geological and tectonic setting

The Zermatt-Saas Zone (ZSZ) and the Theodul Glacier Unit (TGU) are two juxtaposed tectonic units located in the Western Alps (Fig. 1). The ZSZ metaophiolite, which extends from Saas-Fee (Switzerland) in the north to the Aosta Valley (Italy) in the south, represents one of the largest and deepest known portions of exhumed oceanic lithosphere in the world (Bearth 1967; Dal Piaz and Ernst 1978; Reinecke 1998; Bucher et al. 2005; Angiboust et al. 2009). The unit preserves a dismembered ophiolitic sequence made of an association of serpentinites, metagabbros, metabasalts, marbles, schists, calcschists, meta-radiolarites and Mn-bearing metacherts (Bearth 1967; Dal Piaz 1979). The overall lithostratigraphy has slow spreading ridge affinity (Dilek and Furnes 2014) with an internal tectonic structure poorly constrained and highly variable (Angiboust and Agard 2010). The 150–500 m thick ophiolite slices are systematically underlain by a serpentinite body that is 500 m to several km thick (Angiboust and Agard 2010). The rocks have experienced eclogitic high pressure (HP) metamorphism at 23–25 kbar and 530–550 °C (Bearth 1967; Dal Piaz and Ernst 1978; Bucher et al. 2005; Angiboust et al. 2009), with small lenses showing typical ultra-HP mineral relicts corresponding to a pressure of 26–28 kbar (Reinecke 1991; Groppo et al. 2009). The metamorphic ages range between 50 and 40 Ma (Rubatto et al. 1998; Lapen et al. 2003; Dragovic et al. 2020). Retrogression during exhumation locally reflects rock re-equilibration under blueschist facies conditions at 15–16 kbar (Angiboust et al. 2009).

Fig. 1
figure 1

a Geological map of the TGU. A-A’ specifies the location of the cross section (see Online Resource 1: Fig. S1). White circles show sample location. b Inset showing the position of the TGU within the Western Alps

The TGU is located south of Zermatt in Switzerland (Fig. 1) and is a 2 km2 and 100 m thick tectonic slice made of complex interlayering between mafic and felsic rocks interpreted as a volcanoclastic sequence. The main rock types are mafic fels, mafic schist, chloritoid-schist and garnet-schist (Cld-schist and Grt-schist hereafter). The TGU is now folded and fully embedded within the ZSZ (Online resource 1: Fig. S1). The complex garnet textures in the mafic schist reflect a multi-stage metamorphic history of the TGU, which is under debate. Bucher et al. (2019, 2020) interpreted the TGU as a polymetamorphic continental slice based on thermodynamic modelling of garnet in the schists. This interpretation would equal the TGU to the other slivers of continental crust that are found within the ZSZ (e.g. Emilius, Etirol Levaz, Verres, Touir Ponton, Glacier-Rafray, Acque Rosse). Bovay et al. (2021) conducted Lu–Hf dating on garnet from the schists and mafic fels of the TGU and obtained ages ranging from 48 to 50 Ma, which favoured the interpretation of a mono-metamorphic Alpine history. This age range corresponds to the metamorphic ages previously determined for ZSZ (Rubatto et al. 1998; Amato et al. 1999; Dal Piaz et al. 2001; Lapen et al. 2003; Dragovic et al. 2020) and might represent the early stage of Piemont-Ligurian oceanic lithosphere subduction as suggested by Weber et al. (2015) and Bovay et al. (2021). Thermodynamic modelling and Zr-in-rutile thermometry additionally show that TGU mafic fels and metasediments reached a common maximum metamorphic condition at 26.5 ± 1.5 kbar and 580 ± 15 °C (Bovay et al. 2021), which is similar to that reported for the ZSZ (Bearth 1967; Dal Piaz and Ernst 1978; Bucher et al. 2005; Angiboust et al. 2009). In the metasedimentary rocks, an external garnet rim with sharp chemical zoning was assigned to a reheating step during exhumation at ca. 15–17 kbar from 555 to 590 °C (Bovay et al. 2021).

Samples

The three main lithologies of the TGU were investigated: mafic fels, mafic schist and the Cld-schist. The samples contain a rather typical assemblage for blueschist to eclogitic rocks and thus are not identified as metasomatic rocks with limited and/or unusual mineralogy. Veins within the mafic fels are represented by sample Z17TB04A and Z17TB05v (Table 1). Multiple samples were collected over the area to get a spectrum and broad representation of rock-type variability within the tectonic unit (Table 1). These lithologies are described in more details in Bovay et al. (2021) and briefly summarized below.

Table 1 Sample list

Mafic felses outcrop either as large eclogitic bodies or as stretched mafic boudins that are embedded within the surrounding schists (Fig. 2a) and elongated parallel to the main foliation. The thickness of the mafic bodies varies from 1 to 10 m, whereas the mafic boudins have various thickness and elongations ranging from some centimetres to pluri-metres in size. Two end-members mineral associations are observed: (1) garnet-omphacite-glaucophane-rutile-quartz-phengite (sample Z18TB15); (2) garnet-diopside-clinozoisite-zoisite-titanite-graphite. Thermodynamic simulations support the stability of assemblage (1) at HP conditions within the eclogite-facies; zoisite-clinozoisite-rich domains are produced through lawsonite breakdown reaction upon initial exhumation (Bovay et al. 2021). In both mafic fels assemblages, garnet crystals are mostly euhedral having a diameter between 0.2 and 1 cm. Some veinlets were found in zoisite-clinozoisite-rich mafic boudins (sample Z17TB04A, Z17TB05v). In some veins, a symmetrical continuous layer of euhedral garnet crystals can be observed at the contact with the host rock, whereas the inner part of the vein consists of diopside-amphibole-quartz-calcite (Fig. 2b). These minerals do not show any preferred orientation. Amphibole with little quartz often found around diopside is interpreted as a retrogression product.

Fig. 2
figure 2

Field images of TGU lithologies. Mineral abbreviation are from Whitney and Evans (2010). a Mafic boudin wrapped by the foliation within mafic schist. The fine interlayering of the schist has various compositions (see text). b Vein within epidote-rich mafic fels. cd Mafic schist with numerous pseudomorphs after lawsonite having the typical diamond shape. The detail d shows different types of lawsonite relict area, i.e. with and without garnet (see text). e Cld-schist with large flakes of dark chloritoid

The mafic schist (sample Z16TB04, Z16TB11, Z16TB28) is banded and foliated, with layers that range from a few millimetres to a few metres in thickness (Fig. 2a). The common mineral assemblage is garnet-phengite-paragonite-quartz, whereas zoisite, clinozoisite, chloritoid and amphibole complete this mineral association in various proportions. Chloritoid never exceeds 1–2 vol%, whereas garnet, zoisite, clinozoisite and amphibole have highly variable abundance up to 70%. Irregularly distributed lawsonite pseudomorphs (0–35 vol%) occur as milky euhedral diamonds of circa 0.5–2 cm in size, which are stretched along the main foliation (Fig. 2c). The pseudomorphs are filled with paragonite, zoisite and quartz, with rare garnet (Fig. 2d). Three main garnet populations were characterized. (i) Garnet clusters, composed of grains that are up to 200 µm in diameter. They are either stretched along the main foliation or form continuous millimetre thick garnet-rich layers. (ii) Euhedral garnet grains, which size range between 200 and 700 µm in diameter. (iii) Euhedral garnet porphyroblasts with a typical diameter of 2–7 mm.

The Cld-schist shows a pervasive foliation and is primarily composed of quartz, phengite, paragonite, chloritoid, garnet, amphibole and rutile (sample Z16TB24A). Garnet grains are euhedral and equally distributed with an abundance of 20 vol%. The typical garnet diameter varies from 200 to 400 µm with episodic occurrences of atoll garnet.

Analytical methods

Electron probe micro-analysis and chemical mapping

Electron probe micro-analysis (EPMA) of garnet was performed using a JEOL JXA-8200 superprobe at the Institute of Geological Science, University of Bern. Spot analyses were performed using 15 keV accelerating voltage, 20 nA specimen current, and 40 s dwell time (10 s for each background after 20 s on peak). Nine oxide components were obtained using synthetic and natural standards: almandine (SiO2, Al2O3, FeO), albite (Na2O), anorthite (CaO), orthoclase (K2O), forsterite (MgO), ilmenite (TiO2) and tephroite (MnO).

Quantitative compositional maps were generated from X-ray intensity maps using spot analyses acquired in the same area as internal standards. The X-ray maps were measured by WDS with 15 keV accelerating voltage and 100 nA specimen current, with various dwell time and resolutions. Nine elements (Si, Ti, Al, Fe, Mn, Mg, Na, Ca, K) were measured at the specific wavelength in two successive passes (with Na, Ca and K analyzed on the first scan). Compositional maps were processed using XMapTools 3.2.1 (Lanari et al. 2014, 2019). Representative compositions of each garnet and phengite growth zones were obtained by averaging pixels from manually selected areas. The domains were chosen in a way to avoid any mixing between garnet generations and excluding mineral inclusions.

Oxygen isotopes

In situ oxygen isotope analysis of garnet was performed using a CAMECA SIMS 1280HR ion probe at the Institute of Earth Sciences, University of Lausanne. The analyses were carried out following the protocol detailed in Vho et al. (2020b) and using the primary standard UWG2 garnet (δ18O = 5.80‰, Valley et al. 1995). Count rates on the 16O peak were between 1.2 and 2.0 × 109 counts per second. An off-line matrix bias correction based on the grossular and spessartine contents of garnet was applied using the equation given in (Vho et al. 2020c). In every session, the secondary garnet reference material GRS-JH2 (Xgrs = 0.833 ± 0.008; Vho et al. 2020b) was analyzed to check the validity of the matrix bias calibration (Online Resource 2). The δ18O value obtained for the secondary reference material was systematically within the reference values reported in Vho et al. (2020b). The garnet chemical composition of each domain needed for matrix correction was taken from EPMA analyses acquired prior to SIMS analyses. The resulting matrix bias correction on the δ18O of individual analyses ranges typically between 0.1 and 1.5‰, and is mainly or exclusively related to the grossular component. The internal uncertainty on individual oxygen isotopic analyses ranges between 0.14 and 0.30‰ (2σ), whereas the total uncertainty that includes the repeatability of the primary standard and the residuals on the matrix correction curve is 0.31–0.40‰ (2σ). Therefore, the final uncertainty on average values is forced to be no less than 0.3‰ to account for accuracy.

Forward modelling

Fluid-rock interaction during rock metamorphic history was investigated using the computer program PTLOOP (Vho et al. 2020a), which combines equilibrium thermodynamic calculations with oxygen isotopic fractionation modelling. The model performs successive Gibbs energy minimizations computed using Theriak-Domino (De Capitani and Brown 1987; De Capitani and Petrakakis 2010) along a given PT trajectory and calculates at each simulation step the oxygen isotopic fractionation between the mineral predicted to be stable. The model geometry was set to a two-rock stack with the possibility for any free fluid released by the underneath lithology to migrate and interact with the lithology situated above. As discussed in Vho et al (2020a), an additional external free fluid phase can be introduced in the lower lithology at any simulation step.

Oxygen isotopic variation among stable phases, bulk and free fluid was monitored using an internally consistent database for oxygen isotope fractionation (Vho et al. 2019). Thermodynamic modelling was based on the internally consistent thermodynamic dataset of Holland and Powell (1998) and subsequent updates compiled in tc55. The chemical system was restricted to the system SiO2–TiO2–Al2O3–FeO–MnO–MgO–CaO–Na2O–K2O–H2O–O2. The following a–X relations for solid solutions phases were used: calcite-magnesite-dolomite (Holland and Powell 2003); feldspar (Holland and Powell 2003; Baldwin et al. 2005); spinel, biotite (White et al. 2007); epidote, cordierite, talc, staurolite, chlorite, carpholite and garnet (Holland and Powell 1998); clinopyroxene (Green et al. 2007); chloritoid (White et al. 2000); white mica (Coggon and Holland 2002); amphibole (Diener et al. 2007).

Results

Garnet major element chemistry

The diverse TGU rock types contain garnet with a wide range of textures and chemical compositions. Some key garnet crystals from samples used in this study (mafic fels Z16TB32 and Z18TB15, mafic schist Z16TB11 and Cld-schist Z16TB24A) were already investigated for PT modelling and geochronology in Bovay et al. (2021), yet additional grains and samples are considered here. Representative element maps are shown in Figs. 3, 4 and 5. The end-member fractions of garnet for all samples and lithologies are summarized in Online Resource 1 (Fig. S2). Representative garnet compositions are presented in Online Resource 3. Additional EPMA quantitative maps and BSE images are presented in the Online Resource 1 (Fig. S3, S4).

Fig. 3
figure 3

Representative compositional map of garnet from for the mafic fels samples showing garnet end-members (Xgrs − Xalm − Xprp − Xsps), processed with XMapTools (Lanari et al. 2014, 2019). On the Xgrs map, average garnet δ18O values are superimposed with the location of SIMS spot analyses (white circles)

Fig. 4
figure 4

Representative compositional map of garnet from the vein showing garnet end-members (Xgrs − Xalm − Xprp − Xsps). On Xgrs map, average garnet δ18O values are superimposed with SIMS spot analysis location (white circles)

Fig. 5
figure 5

Representative compositional map of garnet from mafic schists (al) and the Cld-schist (mp) showing garnet end-members (Xgrs − Xalm − Xprp − Xsps). On Xgrs map, average garnet δ18O values are superimposed with SIMS spot analysis location (white circles). ad Sample Z16TB11, garnet porphyroblast. eh Sample Z16TB28, intermediate size garnet. il Sample Z16TB04, small garnet. mp Sample Z16TB24A, garnet in the matrix

Garnet cores from the mafic fels are systematically enriched in Mn (Sps14-15, sample Z17TB07 and Z18TB15). In samples Z16TB32, Z17TB05 and Z17TB07 three chemical domains can be recognized (Fig. 3). From core to mantle, both Fe and Mg concentrations (corresponding to Xalm and Xprp in Fig. 3) increase radially from the core towards the rim at the expenses of Ca (Xgrs). A rather homogeneous garnet rim of maximum 150 µm thickness is marked by a decrease in Fe and Mg coupled with an increase in Ca compared to the inner domain. This transition between mantle and rim is discontinuous and can either be sharp or more gradual over a distance of 50 µm. In the core of garnet from mafic fels Z16TB32, zoisite inclusions are highlighted by diffuse Ca-rich halos (Online Resource 1: Fig. S4a, e). Major element zoning in sample Z18TB15 is similar to the previous samples except for Ca, which shows low Ca inner core, a mantle zone enriched in Ca and a low Ca rim. A discontinuous and distinct rim of 200–300 μm enriched in Mg and depleted in Mn, Fe and Ca is visible in part of the crystal (Online Resource 1: Fig. S4f–i).

Garnet crystals located at the contact between vein and metamafic boudin (samples Z17TB04A and Z17TB05v, Fig. 4) have more complex textures. The overall garnet chemical zoning variation from core to rim is comparable to garnet from the host rock (see above). Additionally, the first internal domain is characterized by Mn-rich areas with euhedral shape trailing parallel to the vein direction. This feature is probably reflecting loci of garnet nucleation localization at the time of vein formation (Hollister 1966). The domain D1 (see Fig. 4b) is characterized by a rimward decrease in Mn and Ca, whereas Mg and Fe increase. Two garnet domains have symmetrically overgrown the first domain perpendicular to the vein direction. Compared to the first domain, in the second domain D2, Fe decreases while Ca increases rimward. The last overgrowth corresponding to domain D3 is wider on the host rock side and shows a dramatic increase in Ca coupled with a decrease in Mg, Fe and Mn. The transition between D2 and D3 is rather sharp, whereas D1 and D2 are cross cut by discordant channels having the chemical composition of the garnet from D3 (Fig. 4).

Garnet grains from the mafic schists (sample Z16TB04, Z16TB11 and Z16TB28) can be classified in three distinct types, which correlate with grain-size variation. Type I garnet consists of porphyroblasts having a diameter over 2 mm. These garnet crystals have a continuous bell-shape growth zoning (Fig. 5a–d), where Fe, Ca and Mn contents decrease from core to rim at the expenses of Mg. Some garnet grains display irregularly distributed Ca-rich patchy areas ranging from 10 to 100 µm in size. Type II garnet forms grains with diameter between 200 and 700 µm (Fig. 5e–h). This garnet type shows a complex chemical zoning and has a highly variable texture with decoupling between major elements. Both Ca and Mn show similar zoning patterns with decreasing contents from core to rim. By contrast, Fe and Mg chemical zoning is more irregular and complex. A discordant and irregular rim of maximum 100 µm thickness is separated from the core by a sharp chemical transition over 1–3 µm. Compared to the garnet core, the rim is depleted in Ca and enriched in Mn and Mg. The chemical compositional range of this garnet type is similar to that of the garnet porphyroblasts (Online Resource 1: Fig. S2). Type III garnet consists of small grains (< 300 µm) and is characterized by a sharp (< 1 µm) and discordant core-rim texture (Fig. 5i–l). The xenomorphic core is enriched in Mg, Fe and Mn and depleted in Ca compared to the rim. Unequally distributed irregular Mn-rich zones occur along the contact within the garnet rim (Fig. 5l). The garnet rim has limited zoning with slight decreasing of Ca and Fe coupled with Mg increase. The core chemical composition (Alm64-68Prp23-27Grs13-15Sps1-3) of type III garnet is unique and is never observed for any type I and type II garnet (Online Resource 1: Fig. S2).

Garnet in the Cld-schist (sample Z16TB24a) has a zoning pattern and style that is similar to type III garnet in the mafic schist, with a discordant core-rim texture (Fig. 5m–p). Unlike the mafic schist, in the rim of this garnet Mn is depleted at the contact (compare Fig. 5l and 5p) and then increases outwards. The Mn zoning pattern of the mantle is cross cut by thin radial channel-like textures with rim composition. Additionally, a euhedral external rim of 15 μm thickness that is enriched in Ca and Mn and depleted in Fe and Mg is observed in the the Cld-schist garnet.

Oxygen isotope of garnet

The oxygen isotopic composition of garnet was measured in situ to identify variations in δ18O associated to major element zoning. The results of δ18O analysis are shown in Fig. 6 plotted against Xgrs garnet content as an indicator of chemical zoning. Individual oxygen isotopic measurements are presented in Online Resource 2.

Fig. 6
figure 6

Binary diagram Xgrs vs δ18O for garnet. Garnet δ18O composition was analyzed with SIMS. Garnet composition for matrix correction was taken from the chemical maps (see text). The uncertainty for Xgrs is smaller than the symbol size, whereas the typical uncertainty for δ18O is reported in the legend

Garnet in the different mafic fels samples shows little variation in oxygen isotopic composition, without systematic zoning from core to rim, and consistently low δ18O values (Fig. 3a,e). In the epidote-rich mafic fels sample Z16TB32, values of δ18O are 0.9–2.3‰ in the core, 1.6–2.1‰ in the mantle and 1.2–2.9‰ in the rim. Similar δ18O compositions were measured in the zoisite-rich mafic fels samples Z17TB05 and Z17TB07: values of δ18O are 0.6–1.6‰ in the core, 0.3–1.8‰ in the mantle and 0.8–1.5‰ in the rim. Garnet in omphacite-bearing mafic fels sample Z18TB15 has δ18O values of 0.9–1.2‰ in the core, 1.3–1.8‰ in the mantle, 1.2–1.4‰ in the rim and 1.7–2.1‰ in the external, discontinuous rim. Within each sample, the range of δ18O values for the garnet zones defined by major element zoning is overlapping.

Despite the complex chemical zoning, the δ18O variations measured in the garnet from the vein within the mafic fels samples are small (Fig. 4a). In sample Z17TB04A, the δ18O ranges from 1.4 ± 0.4‰ in the core and mantle to 1.0 ± 0.4‰ in the rim and the discordant areas. Sample Z17TB05v displays comparable δ18O values, where the core value of 1.4 ± 0.3‰ decreases to 0.9 ± 0.3‰ in both directions, i.e. toward the host rock and inner vein side. The cross-cutting garnet domains show lower δ18O values than the surrounding, but similar to the oxygen isotopic composition of the rim, with values ranging from 0.5 to 1.0‰. The extent of δ18O value variation is in accordance with values of garnet from host rock (sample Z17TB05). In both samples, δ18O values are similarly slightly decreasing from core to rim of around 0.4‰.

The mafic schist is the most abundant lithology in the TGU and in every sample (Z16TB04, Z16TB11 and Z16TB28) the small garnet grains of type III show extreme oxygen isotope zoning from a high δ18O core (8.8–11.1‰) to a low δ18O rim (0.2–2.0‰). The step variation in δ18O corresponds to the sharp chemical zoning defined by major elements (Fig. 5i–l). Within the rim itself, the Ca-rich inner part has slightly lower δ18O values than the Mg-rich external part: from 1.7 ± 0.3 to 2.6 ± 0.3‰ in sample Z16TB04 and from 0.8 ± 0.3 to 1.5 ± 0.3‰ in sample Z16TB11. Type II garnet with intermediate size shows moderate δ18O zoning with a slight δ18O increase from core to rim (Fig. 5e), from 1.5 ± 0.3 to 2.3 ± 0.3‰ (Z16TB04), from 0.6 ± 0.3 to 1.0 ± 0.4‰ (Z16TB11) and from 0.2 ± 0.3 to 0.7 ± 0.3‰ (Z16TB28). The chemically discordant rim always shows equivalent δ18O values to the adjacent rim within uncertainty, such as 2.0 ± 0.3‰ (Z16TB04), 1.0 ± 0.3‰ (Z16TB11) and 1.2 ± 0.3‰ (Z16TB28). Type I garnet porphyroblasts in every sample display homogeneous compositions from core to rim (Fig. 5a), with δ18O values of 2.5 ± 0.3‰ (Z16TB04), 1.4 ± 0.4‰ (Z16TB11) and 1.0 ± 0.3 ‰ (Z16TB28). Garnet porphyroblasts and intermediate garnet grains rim have similar oxygen isotopic composition.

The δ18O variation measured in garnet from the Cld-schist is considerable and comparable to what is observed in the mafic schist: a drop of about 7‰ from 10.0 ± 0.4‰ in the core to 3.0 ± 0.4‰ in the rim. Like garnet from the mafic schist (type III), the change in oxygen isotopic composition corresponds to the sharp chemical zoning (Fig. 5m–p).

Discussion

Open system behavior and pervasive fluid flow at high pressure

The PT path followed by the TGU during Alpine subduction was constrained based on equilibrium thermodynamics (Fig. 7) (Weber and Bucher 2015; Bucher et al. 2019; Bovay et al. 2021). For the samples investigated here, the core of small garnet grains from the mafic schist (type III) and Cld-schist was modelled to be stable at 26.5 ± 1.5 kbar and 580 ± 15 °C, whereas the inner part of the rim is stable at 15.5 ± 1.0 kbar and 555 ± 10 °C (Bovay et al. 2021). The sharp chemical boundary between the high-pressure core and the moderate-pressure mantle corresponds to an equally sharp δ18O drop of 7–9 ‰ from garnet core to rim. This sharp zoning is compatible with slow rates of oxygen intergranular diffusion (Vielzeuf et al. 2005; Page et al. 2019) and suggests that the original δ18O composition of garnet core is preserved during exhumation and cooling. It has been demonstrated that, in a closed-system evolution over a temperature range of 150–200 °C, dehydration reactions, fluid loss and mineral fractionation have only minor influence on the mineral δ18O value (< 1‰) and will not influence the δ18O protolith signature (Kohn 1993; Vho et al. 2020a). Consequently, the magnitude of δ18O variation observed in the small garnet grains can only be explained by ingress of external fluids with low δ18O values and in isotopic disequilibrium with the rock, during the TGU first decompression stage. The samples investigated represent the ambient rock type across the unit and were not collected at the contact between the TGU and the ZSZ, where deformation and lithological contrasts are expected to produce even more dramatic fluxes of external fluids (Epstein et al. 2020).

Fig. 7
figure 7

PT–fluid path of the TGU modified from Bovay et al. (2021). Pressure–Temperature path of ZSZ is added as comparison (dashed lines; Angiboust et al., (2009)). The main steps of garnet growth according to Bovay et al. (2021) and corresponding δ18O composition are specified. The δ18O values are reported here as interval, because for each lithology one or multiple sample were analysed. Samples with (?) have a known δ18O composition but unconstrained PT growth conditions. The two potential external fluid sources are plotted: (1) Lawsonite-out reaction in mafic fels sample Z18TB15 simulated using Theriak-Domino (De Capitani and Brown 1987; De Capitani and Petrakakis 2010) assuming water saturation; (2) Atg + Br- > Ol + water range from Johannes (1968). MS mafic schist, MF mafic fels, CldS Cld-schist

Garnet porphyroblasts found in the mafic schist are predicted to be stable during the reheating stage (from 15.5 ± 0.5 kbar and 555 ± 10 °C to 17.0 ± 1.0 kbar and 590 ± 15 °C, Fig. 7) that follows the initial decompression stage (Bovay et al. 2021). Both garnet porphyroblasts (type I) and intermediate garnet with complex textures (type II) display a homogeneous oxygen isotope composition. Their oxygen isotope signature is comparable to that of the rim of small garnet grains. This relative timing indicates that fluid-rock interaction associated to the drop in δ18O occurred prior to the growth of the large garnet grains. Moreover, the small amount (~ 5 vol%) of relic high-δ18O garnet core versus the widespread and high modal abundance (15–20 vol%) of the larger garnet grains suggest that the fluid at the origin of the δ18O drop was pervasive and affected the entire tectonic unit. The mechanism that leads to the complex textures in the intermediate size garnet of type II is unclear (Fig. 5e–h), but it occurred after the ingress of low-δ18O external fluids.

In the mafic fels, garnet formed during the prograde and retrograde evolution of the TGU, and thus the outer rim was already formed during the ingress of external fluids documented in the country rock schists (Fig. 7). Garnet from sample Z16TB32 recorded prograde (490 ± 10 °C and 17.5 ± 0.5 kbar) to peak pressure conditions (580 ± 10 °C and 26.5 ± 0.5 kbar), whereas garnet from sample Z18TB15 recorded both prograde (490 ± 10 °C and 17.5 ± 0.5 kbar) and the reheating stage occurring after the first decompression (590 ± 10 °C and 17.0 ± 0.5 kbar) (Bovay et al. 2021). All the garnet grains analyzed show similar δ18O signature ranging between 1 and 3‰, indicating an unchanged δ18O reactive bulk composition in the mafic fels throughout the entire TGU metamorphic history. The oxygen isotopic composition of the garnet from the mafic fels is equivalent to that of the garnet from the schists, excluding the core of the small garnet grains (type III). Because the prograde-to-peak garnet in the mafic fels already grew from a low δ18O reactive bulk composition, it is impossible to track if this rock type also interacted with the low-δ18O external fluid that affected the metasediments, despite the mafic fels exhibiting a hydrated assemblage with abundant epidote and amphibole.

Veins in metamorphic terrains are commonly interpreted as evidence for fluid flow along fractures, which develop when fluid pressure is close to lithostatic (Ferry 1994). The preserved anhydrous mineral paragenesis (garnet + diopside) observed in veins found within the TGU mafic boudin suggests vein formation during prograde metamorphism at eclogitic facies through devolatilization reactions (Schmidt and Poli 2003). The complex textures observed in garnet from the vein (Fig. 4) point to multiple episodes of garnet growth, associated to (hydro-) fracturing and possibly dissolution–precipitation. Similar processes were documented in other high-pressure terranes (e.g. Angiboust et al. 2012, 2017; Giuntoli et al. 2018a, b; Hyppolito et al. 2019; Rubatto et al. 2020; Vho et al. 2020b). Garnet likely grew in the presence of fluids but did not record significant variation of the reactive bulk δ18O composition. This suggests that the fluids that fractured and partially dissolved the garnet in the vein were in (close) isotopic equilibrium with the rock. Subsequently, the formation of a hydrous retrograde assemblage, requiring more H2O than available from the peak assemblage with amphibole and quartz around diopside, testifies to further rehydration during decompression.

Oxygen isotopic composition of the fluid-rock system and potential fluid sources

The reactive bulk oxygen isotopic composition of each sample can be reconstructed for a given temperature knowing the garnet δ18O value, mineral modal abundances in equilibrium at temperature of garnet growth, and thermodynamic oxygen isotope fractionation factors among coexisting phases. Textural evidence combined with PT data suggest, as discussed above, that the oxygen isotope shift in the mafic schist occurred during the first decompression stage. The TGU unit underwent multiple fluid influxes at the oceanic and subduction stage and thus the original bulk rock δ18O cannot be retrieved from measuring the bulk rock δ18O composition. The original bulk rock δ18O composition was instead calculated for each lithology at peak pressure conditions, before the ingress of external fluids (Online Resource 4), using the software PTLOOP (Vho et al. 2020a) and the bulk rock compositions for major and minor elements (Online Resource 5). The use of bulk rock chemistry for phase equilibrium modelling requires the assumption of minor or insignificant bulk chemical variations, outside water, during fluid-rock interaction (see discussion in Bovay et al. 2021). For the following calculations, a value of 1‰ is taken as average of the range of low δ18O values measured in garnet from all rock types (0–3‰, Fig. 6).

An original bulk δ18O composition of ~ 2‰ is calculated for the mafic fels in equilibrium with the measured garnet δ18O of 1‰. Such low bulk δ18O is below the MORB range (Taylor 1968; Muehlenbachs and Clayton 1972; Alt et al. 1986) but is reported for hydrothermally altered mafic oceanic crust (Cartwright and Barnicoat 1999; Miller et al. 2001; McCaig et al. 2007). Therefore, the relatively low bulk δ18O values of the mafic fels could be ascribed to high-temperature seafloor alteration prior to Alpine metamorphism. A similar process affecting the schists would have equally decreased their original δ18O of a few ‰ in the oceanic stage. For the mafic schist and Cld-schist, the calculated bulk δ18O is ~ 11‰ based on the δ18O of 9.5‰ of the garnet core. The bulk rock δ18O composition of the schists corresponds to the lower side of the wide δ18O range usually observed for metasediments and more specifically pelites (Rumble et al. 1982; Bebout and Barton 1989; Cartwright and Barnicoat 1999). On the other hand, a bulk rock with δ18O of ~ 3‰ would be in equilibrium with the low δ18O garnet rim in the schists (Online Resource 1: Fig. S7, S8).

Similarly, the δ18O value of the fluid in isotopic equilibrium with garnet at the time of (or following) fluid–rock interaction was calculated (Online Resource 4), where fluid is considered as pure H2O. Considering a temperature range of 580–555 °C, pure water with δ18O of ~ 3.5‰ would be in equilibrium with the ~ 1‰ garnet rim from the schists (and coincidentally with the garnet from the mafic fels). This allows us to conclude that the decrease in δ18O in the schists was due to an external fluid with a δ18O of ~ 3.5‰ or lower, as any partial equilibration with a fluid with lower δ18O could produce the same garnet and bulk rock signature. The calculated values for the reactive bulk rock δ18O composition and fluid in equilibrium with the garnet rim are the basis for reconstructing the potential fluid source and the fluid dynamics (see below).

Fluid is commonly thought to be expelled during prograde metamorphism in subduction zones (Schmidt and Poli 2013), and additional fluids can also be generated via retrograde dehydration reactions (Angiboust and Agard 2010). Slab material being exhumed towards the surface during retrogression, such as the TGU, is expected to interact with external fluid either expelled from structurally lower layers within the slab or from deeper down in the subduction system. Prior to subduction, TGU was a volcanoclastic sequence (Bucher et al. 2020; Bovay et al. 2021) and was stratigraphically the uppermost layer atop the magma-poor Piemont-Ligurian oceanic crust, now identified as the ZSZ metaophiolite. Considering the rock types involved and the geological relationships, the fluid that interacted with TGU rocks during early exhumation could be provided by two main mineral dehydration reactions (Fig. 7): (1) lawsonite-out within the TGU or ZSZ mafic rocks; (2) antigorite-out within the underlying hydrated ultramafic rocks (serpentinites) from the ZSZ. For these two scenarios, the reactions involved and the isotopic signature of the fluids are considered below.

(1a) Lawsonite-out dehydration reaction, in TGU. Lawsonite is a common metamorphic mineral at HP–LT conditions and contains 11.5 wt% H2O in its structural formula, which makes it is one of the major water carriers in subduction zones (Schmidt and Poli 1998; Clarke et al. 2006). Preserved lawsonite is seldom observed in metaophiolites even though experimental work and modelling suggest that lawsonite eclogite should be an abundant rock type in oceanic subduction zones (Poli and Schmidt 1995; Clarke et al. 2006; Angiboust and Agard 2010). Instead, lawsonite is often retrogressed in epidote-rich aggregates during decompression and exhumation (Ernst and Dal Piaz 1978; Ballèvre et al. 2003; Angiboust and Agard 2010). Typical diamond shape pseudomorphs after lawsonite are commonly found in the TGU schists, but never in the mafic fels. However, this rock type contains various amounts (up to 70 vol% in sample Z17TB05) of epidote group minerals (Fig. 2b). Assuming water saturation and using the bulk rock chemistry presented in Online Resource 5, a thermodynamic forward model of mafic fels at TGU maximum metamorphic conditions (580 °C, 26.5 kbar; Bovay et al. 2021) predicts lawsonite modal abundances ranging between 9 vol% (sample Z18TB15) and 40 vol% (sample Z17TB05). During exhumation of the TGU, lawsonite breakdown reaction in the mafic fels (low δ18O bulk) will release a H2O fluid having the required δ18O composition to change the bulk δ18O in the TGU schists (see above).

(1b) Lawsonite-out dehydration reaction, in ZSZ. Lawsonite relicts have been described in the ZSZ metabasalts (Angiboust and Agard 2010) that have bulk rock δ18O values from 4.4 to 9.2‰ (Barnicoat and Cartwright 1995; Cartwright and Barnicoat 1999). Such bulk rock compositions would produce fluids with even higher δ18O values. As calculated above, the entering fluids in the TGU metasediments must have had a δ18O value of 3.5‰ or lower. Therefore, the fluid released by lawsonite breakdown in the metabasic rocks from the ZSZ during exhumation (Angiboust and Agard 2010) have a δ18O values that is too high to be the fluid entering the TGU schists. Metabasalts with lower δ18O have never been reported for the ZSZ, whereas ZSZ serpentinites yield δ18O values of 1–3‰ (Cartwright and Barnicoat 1999). Additionally, quantitative fluid modeling excludes the lawsonite-out reaction as source of the necessary amount of external fluid for scenario 1a and 1b (see below).

(2) Antigorite-out dehydration reaction. Subduction of serpentinized peridotite is the most efficient mechanism to carry water at depth (Scambelluri et al. 1995) and the most important dehydration reaction occurring in serpentinite is related to antigorite breakdown (Ulmer and Trommsdorff 1995; Schmidt and Poli 1998). Antigorite dehydration occurs through two subsequent dehydration reactions: (a) antigorite + brucite→olivine + chlorite + water (Johannes 1968); (b) antigorite→olivine + enstatite + water (Ulmer and Trommsdorff 1995; Padrón-Navarta et al. 2013). Reaction (a) occurs within the temperature range ~ 500–550 °C, which overlaps with the PT history of the TGU (Fig. 7), whereas reaction (b) corresponds to the full antigorite breakdown that is predicted to occur at higher temperature (and pressure) in the downgoing slab. As a slow spreading setting, the Piemont-Ligurian oceanic lithosphere was harzburgitic and was largely hydrated at the ocean floor (Ernst and Dal Piaz 1978; Li et al. 2004; Dilek and Furnes 2011). Kempf et al. (2020) showed that partly hydrated harzburgite can contain from ~ 5 to 10 wt% brucite and reaction (a) would release between 3.4 and 7.2 wt% of water at 25 kbar. Literature data show that ultramafic rocks in the Alps have a wide range in δ18O, from 1 to 9‰ (Früh-Green et al. 1990, 2001; Cartwright and Barnicoat 1999; Miller et al. 2001). Serpentine-water δ18O equilibrium calculations, using the fractionation factors of Vho et al. (2019), show that at 500–550 °C serpentine with a δ18O of 2‰ is in equilibrium with water having a δ18O value of 3.7–3.9‰. This value is comparable to the maximum δ18O of the fluid entering the TGU metasediments (see above). The oxygen isotope composition of serpentinite adjacent to the TGU is between 4.1 and 5.2‰ (Lucile Seydoux unpublished data), serpentinites from Mellichen and Valtournanche, which are relatively near to the TGU, have bulk δ18O values between 4.3–5.8‰ and 1.2–3.2‰, respectively (Cartwright and Barnicoat 1999). Serpentinites comparable to those outcropping in Valtournanche could thus represent a suitable source for fluid infiltrating the TGU rocks.

Quantitative fluid modelling

Modelling using the software PTLOOP (Vho et al. 2020a) was conducted to evaluate the validity of the two scenarios introduced above and results are shown in Table 2 and Online Resource 1 (Fig. S5, S6, S7, S8). A decompression PT trajectory from 580 °C and 26.5 kbar to 555 °C and 15.5 kbar (Bovay et al. 2021) including 15 successive simulation steps was used. For each simulation step, the mineral paragenesis, δ18O composition of stable phases, mass and δ18O composition of excess fluid for the source rock are calculated. At every step of the simulation, the totality of free fluid predicted to be stable in the source rock is transferred to the overlying metasediment. After fluid influx, the mineral paragenesis, δ18O composition of stable phases, mass and δ18O composition of excess fluid for the metasediment are calculated. Finally, the excess fluid from the metasediment is removed from the system before the next simulation step begins. Prior to run the simulation, the water content required for the observed paragenesis at maximum metamorphic condition was investigated. Pressure–temperature phase diagrams were calculated with associated pristine or fractionated bulk rock chemistry (Online Resource 5), using the Theriak Domino software (De Capitani and Brown 1987; De Capitani and Petrakakis 2010). In every case, the best result was achieved with water saturation and this assumption was set as an input parameter in PTLOOP (Online Resource 1: Fig. S9a, d). The pervasive δ18O resetting in the schists attests for fluid-rock interaction at the scale of the tectonic unit. The modelling is performed considering solely the mafic schist among the metasedimentary rock, because this is the main metasedimentary type observed within the TGU. For the simulation, the bulk rock chemistry of the mafic schist (Z16TB11) was adapted by fractionating 5 vol% of garnet core (estimate based on petrographic observations), considering that this garnet relict did not react at the peak stage (Online Resource 5).

Table 2 Input and output values from PTLOOP simulations to evaluate the fluid origin and reconstruct the time-integrated fluid flux

The first scenario considers a rock column (height × 1 [m] × 1 [m]) with 100 m of mafic schist on top of 10 m of mafic fels. The thickness ratio between the mafic fels and the mafic schist was estimated based on field observations, where the lithologies have a subhorizontal stacking (Weber and Bucher 2015). The bulk chemical composition of mafic fels Z18TB15 was chosen because it is the most common type of mafic fels found in the TGU. Results (Table 2) show that water expelled from the mafic fels and infiltrating the mafic schist has a δ18O signature of 3.5‰ for a total integrated mass of 161 kg of H2O. This corresponds to a fluid-rock mass ratio of ~ 5.8 × 10–4 assuming a rock density of 2800 kg/m3. Infiltration of this amount of fluid into the mafic schist produces only a minimal shift in garnet δ18O composition from 9.5 to 9.3‰, which is not enough to match the observed shift. To achieve the garnet oxygen isotopic composition of ~ 1‰ (thus − 8.5‰ shift), the mafic schist/mafic fels volume ratio has to be set to an unrealistic value of 1:1500 (Table 2), which is not in line with the proportion of each rock type observed in either the TGU, or in the ZSZ. Therefore, the mafic fels alone cannot produce enough fluid to explain the observed δ18O shift.

The second scenario considers external fluid coming from the serpentinite (a rock column of 100 [m] × 1 [m] × 1 [m], row 3 in Table 2) and infiltrating directly the metasediment. This procedure allows an exact evaluation of the quantity of external fluid needed to shift garnet δ18O composition from ~ 9.5 to 1‰, and the bulk rock accordingly. The external fluid δ18O composition was set at 3.5‰, which is in isotopic equilibrium with the mafic fels and with the range of δ18O measured in the Valtournanche serpentinites (Cartwright and Barnicoat 1999). The ingress of serpentinite-derived fluid was spread over several successive batches to simulate a limited continuous fluid flow. To produce the observed dramatic shift of -8.5 ‰ in garnet and bulk rock, 1.1 × 106 kg of H2O from the serpentinite is required to enter the system (Table 2). This corresponds to a fluid-rock mass ratio of ~ 3.9 assuming a host rock density of 2800 kg/m3. Given that serpentinite releases from 3.4 to 7.2 wt% water through the reaction antigorite + brucite→olivine + chlorite + water (Kempf et al. 2020), a 5.6–11.8 km thick layer of hydrated ultramafic crust (density of 2750 kg/m3) is required below the metasediments to produce the necessary amount of fluid. A serpentinite layer of such thickness would be a maximum estimate. The TGU is volumetrically small compared to the underlying oceanic lithosphere exposed in the ZSZ and serpentinite-derived fluids could be sourced from a much larger area, reducing the thickness of the serpentinite layer needed to provide a sufficient amount of fluid. Hydration of the oceanic lithosphere commonly occurs in the first two kilometres, but partial hydration could reach greater depth of ~ 10 km, particularly in slow-spreading ridges (Schmidt and Poli 1998; Baumgartner and Valley 2001; Früh-Green et al. 2001; Angiboust and Agard 2010; Dilek and Furnes 2011). Therefore, the serpentinite are identified as a suitable source for the fluids in terms of isotopic signature and availability.

Pervasive flow, time-integrated fluid flux, rock permeability

Channelized fluid flow in one of more high permeability conduits, such as veins, fractures and lithological contacts, is regarded as the representative mechanism for intra-slab fluid flow (Miller et al. 2001; Scambelluri and Philippot 2001; Hermann et al. 2006; Zack and John 2007; Ague 2011; Spandler et al. 2011). However, in some parts of subducting slabs, important fluid release could be pervasive, namely fluid migration around individual mineral grains through an interconnected porosity (Ague 2014). This latter behaviour is favoured in units with monotonous rheology and permeability (Oliver 1996). In the TGU schists, the decrease in δ18O from garnet cores at 9–11% to values that converge at 0–3‰ (Fig. 6) testify for pervasive fluid flow across a 2000 × 100 m tectonic unit during HP metamorphism. This is unlike the scenario proposed by studies that describe small-scale channelized fluid flow in networks of veins and shear zones, whereas the country rock shows a “closed-system” behaviour (Zack and John 2007; John et al. 2012; Jaeckel et al. 2018; Epstein et al. 2021). Large fluid fluxes over km scales documented in other localities of the Western Alps metaophiolites have generally been associated to shear zones (e.g. Angiboust et al. 2015; Jaeckel et al. 2018), whereas the TGU is a coherent package of schist and fels.

The time-integrated fluid flux corresponds to the total amount of external fluid influx that passes across an area of interest during a given time interval (Ague 2003). Fluid channelization processes strongly impact time-integrated fluid fluxes (Fig. 8). It has been proposed by Zack and John (2007) that at subarc depth, the time-integrated fluid flux of the H2O escaping the subducting oceanic crust regarded as pervasive fluid flow is ~ 3.4 × 104 cm3/cm2. However, quantifications of time-integrated fluid flux in subduction zone where pervasive flow is proposed remain scarce. This is because open-system behaviour can only be traced in the rare cases where the zones are preserved and exposed to the surface where they become accessible to investigation (Zack and John 2007); the TGU offers this rare opportunity.

Fig. 8
figure 8

Time-integrated fluxes in different geological environment (modified from Zack and John 2007). Circles represent punctual value, whereas rectangles represent range of analyses. The boundary between open and closed system is from Oliver (1996). Fluid fluxes for oceanic lithosphere a eclogitic body is from Zack and John (2007), b pervasive H2O flow at subarc depth is deduced by Zack and John (2007) using data from Schmidt and Poli (1998), and c channelized H2O flow at subarc depth is from Zack and John (2007), where flux intensity is inversely proportional to the vein density. Fluid fluxes for: upper crust is from Ferry and Gerdes (1998), for contact metamorphism is from Ague (2003), accretionary prism is from Ague (2003), eclogite fracture zone is from Philippot and Rumble (2000), and Sesia Zone: white rectangle is from Konrad-Schmolke et al. (2011), and black rectangle is from Vho et al. (2020b)

In the TGU metasedimentary rocks, the time-integrated fluid flux is the total amount of fluid involved in the fluid-rock interaction that leads to the dramatic δ18O drop throughout the unit. It can be deduced from the total amount of fluid expelled from the serpentinite that interacted with the TGU (water density of ~ 103 kg/m3) and corresponds to 1.1 × 105 cm3/cm2 (Table 2). This value represents a conservative estimate, because a fraction of the fluid could pass through the system without reacting with the rock (Zack and John 2007). This time-integrated flux is independent of deformation, which may speed or slow down flow rates at any given time, but does not affected the time-integrated signal. This result is above the time-integrated fluid flux threshold value of 1 × 104 cm3/cm2 defining open-system behaviour (Ferry and Dipple 1991; Dipple and Ferry 1992; Oliver 1996), beyond which a major tracer, like oxygen isotopes, is affected by infiltrating aqueous reactive fluid (Fig. 8).

The time-integrated fluid flux in the TGU can be compared to values proposed for the high pressure continental unit of the Sesia Zone (Western Alps), which underwent a PT evolution similar to the TGU (Compagnoni 1977; Gosso 1977) and was affected by pervasive rehydration during exhumation (Konrad-Schmolke et al. 2011). Time-integrated fluid flux calculated for the Sesia Zone by Konrad-Schmolke et al. (2011) are 1.5–7.5 × 103 cm3/cm2. Such values cannot be classified as open-system behavior according to the definition of Oliver (1996) and are comparable to the time-integrated fluid flux of dehydrated eclogites (5 × 10–5 × 103 cm3/cm2) situated away from channelization drainage (Zack and John 2007). By contrast, Vho et al. (2020b) calculated a time-integrated fluid flux for the Sesia Zone (1.7–2.2 × 105 cm3/cm2) that is above the threshold value for open-system behavior. Likewise, systems such as fracture zones in eclogite (Philippot and Rumble 2000), upper crustal fracture zones (Ferry and Gerdes 1998), contact metamorphic aureoles and accretionary wedges (Ague 2003) typically have high time-integrated fluid fluxes of > 104 cm3/cm2. The above estimates suggest a diversity of fluid fluxes in the slab-mantle transition zone (Konrad-Schmolke et al. 2011).

Fluid migration must operate along distances of several kilometres in the subduction channel (John et al. 2012) and needs highly interconnected porosity to act efficiently (Ague 2003). Parameters taken in account to calculate rock permeability are mostly integrated with respect to time, because fluid-rock interaction take place throughout a discrete time window and must have been sufficient to accommodate fluid flow recorded by the exhumed metamorphic rocks (Manning and Ingebritsen 1999; Ague 2003). Rock permeability can be calculated as follow (Ingebritsen and Manning 1999; Manning and Ingebritsen 1999):

$$k=\frac{Q\hspace{0.2cm} \cdot\hspace{0.2cm}\mu }{\Delta t\hspace{0.1cm} \cdot\hspace{0.1cm} g \hspace{0.1cm} \cdot\hspace{0.1cm}({\rho }_{\mathrm{rock}}-{\rho }_{\mathrm{fluid}})}$$
(1)

where Q is the time-integrated fluid flux (m3/m2) (Table 2), µ is the fluid viscosity (~ 1 × 10–4 Pas at 550 °C for water), Δt is the duration of fluid flow in seconds, g is the gravitational force equivalent (9.81 ms−2), ρrock the rock density (2800 kgm−3 for metasediments), ρfluid the fluid density (997 kgm−3 for pure water). In the deep crust, the maximum fluid pressure are approximatively lithostatic and for vertical upwards fluid flow path the driving-force gradient is approximated by · (ρrock − ρfluid), which is the difference of the lithostatic pressure and hydrostatic pressure (Manning and Ingebritsen 1999; Ague 2003). In the case of the TGU, the duration of fluid–rock interaction is not well constrained, but most geochronological data for HP metamorphism in the ZSZ are within a few million years (e.g. Dragovic et al. 2020). Therefore, a time parameter Δt between 1 and 10 Ma was considered: this leads to permeability variations of a factor 10, where longer duration will decrease the permeability (Fig. 9).

Fig. 9
figure 9

Evolution of permeability with depth compared to results for this study. Empirical formula (bold line) from Ingebritsen and Manning (1999). Grey-shaded squares are different rock permeabilities from experiments on ultramafic rocks (Ganzhorn et al. 2019). Impermeable rock domain (hatched area) from Ganzhorn et al. (2019)

The calculated permeability of the TGU metasediment according to Eq. (1) is 2 × 10–20 m2 for 10 Ma of external fluid influx and 2 × 10–19 m2 for 1 Ma (Fig. 9). An alternative way to calculate permeability is using the empirical formula that relates permeability with depth: log k = -14–3.2 log(z), where z is the depth in kilometres (Manning and Ingebritsen 1999; Ingebritsen and Manning 2010). The permeability calculated with formula (1) for the TGU matches the empirical value at 60 km depth for the 10 Ma scenario, or at 30 km depth for the 1 Ma case (Fig. 9). Simulations indicate that fluid-rock interaction in the TGU occurred between 26.5 and 15.5 kbar. Using a typical continental crust rock density of 2800 kg/m3 this corresponds to a depth range of 95 and 55 km depth and favours a scenario lasting at least 8 Ma (Fig. 9). Because porosity and permeability evolve with deformation (Ague 2014), these remain indicative values only. It has to be noted that the investigated samples were not collected within shear zones or at the localized tectonic contact between the TGU and the ZSZ. It is also noteworthy that our result is comparable to the permeability that was experimentally determined for blueschists and deformed antigorite serpentinite, whereas Atg-serpentinite and chlorite schist show lower permeabilities (Kawano et al. 2011; Ganzhorn et al. 2019).

Fluid flow initiation

Volume variations caused by mineral reactions in the different lithologies of the TGU were investigated to assess fluid flow initiation from the serpentinite of the ZSZ towards the schists of the TGU. Rock volume variation in the TGU mafic fels sample Z18TB15 and mafic schist Z16TB11 was modelled using Theriak-Domino (De Capitani and Brown 1987; De Capitani and Petrakakis 2010). Water saturation is assumed at maximum metamorphic condition and free fluid is not allowed to exit the system. As no evidence of fluid-rock interaction were observed along the prograde path (Fig. 7), simulations were consequently focused on the early stage of decompression from 580 °C and 26.5 kbar to 555 °C and 15.5 kbar. We defined the cumulative rock volume variation as:

$${V}_{\mathrm{cumulative},i}=\left(\frac{{V}_{i}-{V}_{i-1}}{{V}_{0}}+{V}_{\mathrm{cumulative},i-1}\right)\times 100\hspace{0.2cm}, \hspace{0.2cm} where \hspace{0.2cm} {V}_{cumulative,0}=0 \hspace{0.2cm} and \hspace{0.2cm} i\in {\mathbb{N}}$$
(2)

where i corresponds to a specific iteration step along the retrograde path and V0 is the total volume of solids at maximum pressure taken as a reference. Results plotted against pressure show an overall absolute volume increase of ~ 4% along the exhumation path for both lithologies (Fig. 10). Localized volume variations of various magnitude are associated to mineral reactions and occur at different iterative steps (i.e. P and T) in each lithologies (Fig. 10a–d). Interestingly, negative volume variation in both rock types is associated to lawsonite dehydration reaction (Fig. 10, Online Resource 1: Fig. S9). Following the exhumation path, this reaction is predicted to occur first in the mafic fels at ~ 20.0 kbar and later on in the mafic schist at ~ 19.5 kbar, where the rock volume variation is ~ 1% in the mafic fels and ~ 2% in the mafic schist, respectively. Mafic fels from the TGU are not in contact with the ZSZ (Fig. 1 and Online Resource Fig. S1) and show no evidence for external fluid input. Consequently, volume variations in the mafic fels will not perturb the overall fluid flow dynamic during fluid transfer between the ZSZ serpentinite and the TGU schist. Moreover, field observations show that the mafic boudins within the schists are not interconnected. These observations make the possibility of fluid infiltration along structural features, such as lithological contacts or boudins necking unlikely (Fig. 10c).

Fig. 10
figure 10

Cumulative volume variation for selected TGU rock types plotted along the TGU retrograde metamorphic path. Main volume variations are related to hydrous phases breakdown (see details in text). Sketches on the right ad are a representation of the system at the point of hydrous phases breakdown during TGU exhumation. Sketches are not to scale and represent the volume (small black arrows) and porosity (textures) variation of the different lithologies from the TGU and the effect on fluid-rock interaction. Fluid flow is represented by blue lines with arrows indicating the direction of flow. MS mafic schist, MF mafic fels, UM ultramafic

The TGU is mainly composed of mafic schist, and its significant absolute volume variation is assumed to have a major impact on the fluid flow dynamic within the tectonic unit (Fig. 10d). Negative volume changes associated to mineral reactions are known as an efficient mechanism to develop transient water-filled micro-porosity (Rumble and Spear 1983; Oliver 1996; Balashov and Yardley 1998; Ague 2003). Such an event will create short-lived excess permeability needed to let external fluid coming-in and initiate fluid flow. Therefore, lawsonite breakdown in the mafic fels is not a major fluid source in the δ18O budget due to its negligible rock volume, but lawsonite breakdown in the schist is rather a key reaction for the onset of fluid flow from the serpentinites towards the TGU metasediments.

Conclusion

This study identifies the large amounts of low δ18O fluid generated from antigorite + brucite breakdown reaction in the ZSZ serpentinites as the major origin for the external fluid that interacted pervasively with the schists of TGU during exhumation between 26.5 and 15.5 kbar in the subduction channel. Mafic fels from the TGU have no δ18O evidence for an external fluid input and simulations show that they play only a minor role in the δ18O budget of the TGU. Serpentinites from the ZSZ would have the right oxygen isotopic composition to produce the low δ18O fluid that infiltrated the TGU rocks. In the mafic fels, the low bulk δ18O signature of all garnet generations, including those grown during prograde metamorphism, suggests a pre-metamorphic lowering of the bulk rock δ18O, probably through seafloor alteration.

The time-integrated fluid flux calculated for fluid-rock interaction at HP in the TGU falls in the range of open-system behavior. In HP environments, such large fluid fluxes are usually attributed to channelized fluid flow or fracture zone and seldom for pervasive fluid flow as that documented in the TGU. The calculated permeability for the TGU schist during fluid–rock interaction at depth is comparable to values from experimental and empirical studies. Transient rock volume variation associated to lawsonite-out reaction in the schist of the TGU are interpreted to have initiated fluid transfer between the TGU and the surrounding ultramafic rocks from the ZSZ.