Articles | Volume 11, issue 2
Research article
10 Mar 2020
Research article |  | 10 Mar 2020

Tracing fluid transfers in subduction zones: an integrated thermodynamic and δ18O fractionation modelling approach

Alice Vho, Pierre Lanari, Daniela Rubatto, and Jörg Hermann

Oxygen isotope geochemistry is a powerful tool for investigating rocks that interacted with fluids, to assess fluid sources and quantify the conditions of fluid–rock interaction. We present an integrated modelling approach and the computer program PTLoop that combine thermodynamic and oxygen isotope fractionation modelling for multi-rock open systems. The strategy involves a robust petrological model performing on-the-fly Gibbs energy minimizations coupled to an oxygen fractionation model for a given chemical and isotopic bulk rock composition; both models are based on internally consistent databases. This approach is applied to subduction zone metamorphism to predict the possible range of δ18O values for stable phases and aqueous fluids at various pressure (P) and temperature (T) conditions in the subducting slab. The modelled system is composed of a mafic oceanic crust with a sedimentary cover of known initial chemical composition and bulk δ18O. The evolution of mineral assemblages and δ18O values of each phase is calculated along a defined PT path for two typical compositions of basalts and sediments. In a closed system, the dehydration reactions, fluid loss and mineral fractionation produce minor to negligible variations (i.e. within 1 ‰) in the bulk δ18O values of the rocks, which are likely to remain representative of the protolith composition. In an open system, fluid–rock interaction may occur (1) in the metasediment, as a consequence of infiltration of the fluid liberated by dehydration reactions occurring in the metamorphosed mafic oceanic crust, and (2) in the metabasalt, as a consequence of infiltration of an external fluid originated by dehydration of underlying serpentinites. In each rock type, the interaction with external fluids may lead to shifts in δ18O up to 1 order of magnitude larger than those calculated for closed systems. Such variations can be detected by analysing in situ oxygen isotopes in key metamorphic minerals such as garnet, white mica and quartz. The simulations show that when the water released by the slab infiltrates the forearc mantle wedge, it can cause extensive serpentinization within fractions of 1 Myr and significant oxygen isotope variation at the interface. The approach presented here opens new perspectives for tracking fluid pathways in subduction zones, to distinguish porous from channelled fluid flows, and to determine the PT conditions and the extent of fluid–rock interaction.

1 Introduction

The subducting oceanic slab is composed of a sequence of rock types corresponding to chemical systems that undergo continuous and discontinuous reactions in response to pressure (P) and temperature (T) changes. Through its metamorphic history, the hydrated oceanic lithosphere undergoes extensive dehydration by the breakdown of low-temperature, volatile-rich minerals (e.g. Baumgartner and Valley, 2001; Baxter and Caddick, 2013; Hacker, 2008; Manning, 2004; Page et al., 2013; Poli and Schmidt, 2002). The expelled aqueous fluid migrates through the slab towards the slab–mantle interface, and it may continue rising to the mantle wedge playing a major role in triggering mass transfer and melting (Barnicoat and Cartwright, 1995; Bebout and Penniston-Dorland, 2016). Evidence for fluid circulation in subducted rocks has been extensively observed in exhumed high-pressure/ultra-high-pressure (HP/UHP) terrains (e.g. Zack and John, 2007; Baxter and Caddick, 2013; Martin et al., 2014; Rubatto and Angiboust, 2015; Engi et al., 2018), but a direct link to the primary source production is often missing and the main source remains matter of debate. The characterization of fluid pathways in subduction zones has been addressed by using a variety of methods (i.e. seismicity, thermodynamic modelling, fluid inclusions, HP veins, trace element and stable isotope studies on metamorphic minerals) (e.g. Baxter and Caddick, 2013; Hacker, 2008; Hernández-Uribe and Palin, 2019; Scambelluri and Philippot, 2001; Spandler and Hermann, 2005). In particular, oxygen isotope composition of metamorphic minerals from exhumed HP rocks sheds light on the nature of the fluid reacting with those systems during metamorphism. Thus, oxygen isotope studies of HP rocks have the potential to make important contributions to the investigation of fluid sources and pathways in subduction zones (e.g. O'Neil and Taylor, 1967; Muehlenbachs and Clayton, 1972; Baumgartner and Valley, 2001; Page et al., 2013; Martin et al., 2014; White and Klein, 2014; Hoefs, 2015; Rubatto and Angiboust, 2015).

Figure 1Schematic geometry of the subduction models discussed in the text. The rock column is composed of two rock types (Rock 1 and Rock 2) that can be infiltrated by an external fluid deriving from a third layer located beneath them. (a) Example columns used in the calculation along the PT path shown in Fig. 2 to produce the results presented in Figs. 3–7. See text for details. (b) Schematic representation of the three interaction cases discussed in the text. No-interaction case (NI): the fluid released by the metabasalt does not interact with the metasediment. The fluid leaving the system is a mixture of metabasalt-derived and metasediment-derived fluids. Partial-interaction case (PI): 50 % of the metabasalt-derived fluid does not interact with the metasediment, and 50 % of it equilibrates with the metasediment. The final fluid released by the system is the mixture between the unmodified metabasalt-derived fluid and the fluid deriving from the metasediment after it equilibrates with 50 % of the metabasalt-derived fluid. High-interaction case (HI): all the fluid released by the metabasalt equilibrates with the metasediment. The fluid leaving the system exits the metasediment. (c) Possible scenario at the base of the column. As a consequence of serpentine breakdown, serpentinite-derived fluids may infiltrate the metabasalt, exchange with it and affect the fluid infiltrating the metasediment.


The modelling of oxygen isotopic fractionation has been traditionally addressed as an equilibrium calculation between individual mineral couples. An alternative approach follows what has been extensively adopted in the last decades for thermodynamic modelling (see reviews by Lanari and Duesterhoeft, 2019; Powell and Holland, 2008; Spear et al., 2017) and considers an evolving mineral assemblage. A pioneer model proposed by Kohn (1993) can be applied to single and closed chemical systems, i.e. for which no infiltration of external fluids in isotopic disequilibrium is allowed. Such an approach can simulate how the oxygen isotopic composition changes with P and T, but it remains too simple for subduction zone settings, where significant fluid exchange occurs between different lithologies within the subducting slab. Baumgartner and Valley (2001) proposed a model for stable isotope fluid–rock exchange based on continuum mechanics, where infiltration profiles can be calculated, but no information is provided about the different components (minerals) of the rock, as it is regarded as a continuum. The fluid  rock (FR) ratios obtained with this strategy do not correspond to the physical amount of fluid but rather represent a measurement of exchange progress.

We present a new approach that combines equilibrium thermodynamics and oxygen isotope fractionation modelling applied to multi-rock systems. This modelling technique takes advantage of the increased capability of forward modelling of complex systems achieved in the last 2 decades (Lanari and Duesterhoeft, 2019). A MATLAB©-based modelling program PTLoop has been developed to calculate oxygen isotope fractionation between stable phases from the results of Gibbs energy minimization performed by Theriak-Domino (de Capitani and Brown, 1987; de Capitani and Petrakakis, 2010) along any fixed PT trajectory. The oxygen isotope variation in each mineral within the evolving assemblage is tracked using an extensive and internally consistent database for oxygen isotope fractionation (Vho et al., 2020). A graphical user interface (GUI) provides the representation of the results. The capabilities of this software solution will be discussed in detail with an example that focuses on the characterization of (1) the effect of the dehydration reactions on the bulk δ18O of a rock, (2) the effect of the influx into a subducting rock of an external fluid of distinct isotopic composition and (3) the final amount and isotopic composition of the fluid leaving the multi-rock system, e.g. infiltrating an upper unit or the mantle wedge. Petrological implications of relevant computational results are also discussed.

2 Modelling

2.1 Model geometry

The subducting oceanic lithosphere is typically composed of a section of igneous oceanic crust with its sedimentary cover (mostly < 1 km) above and an ultramafic lithospheric mantle section beneath. The geometry of the model is illustrated in Fig. 1. The target column represents a simplified section of the upper part of such oceanic lithosphere. It is composed of a layer of basaltic composition (Rock 1) overlaid by a layer of sediments (Rock 2, see below for details). Two different rock columns are considered: (1) a relatively water-rich system with altered basalts and terrigenous sediments and (2) a relatively water-poor system with unaltered basalts (of mid-ocean ridge basalt (MORB) composition) and carbonate sediments. The column has a fixed section of 1 m2, while the thickness of each rock unit can be set by the user. The model is conservative with respect to the mass, while the volume of each rock type changes according to fluid loss and density variation along the PT path. The PT structure of subduction zones depends on numerous variables, including the age of the incoming lithosphere and the amount of previously subducted lithosphere (e.g. Peacock, 1990). In this study, the calculation was performed following the subduction geotherm from Gerya et al. (2002) (Fig. 2) over a pressure range of 1.3–2.6 GPa, corresponding to a depth of ∼45 to ∼85 km to encompass the conditions of interest for the investigated processes. The modelled temperatures range from a minimum of 350 C to a maximum of 700 C. The lower limit was chosen to take into account the large uncertainties in the thermodynamic databases for many common low-T metamorphic minerals that lead to unsatisfactory models for phase equilibria and mineral parageneses for low-grade rocks (Frey et al., 1991; Vidal et al., 2016). The upper limit is fixed by the wet solidus of metasediments – the present model strictly applies to subsolidus conditions.

Figure 2PT diagram showing typical oceanic subduction geotherms and the intermediate geotherm used in the calculation (red line, Gerya et al., 2002). The red spots represent the modelling steps. The average D80 geotherm from Syracuse et al. (2010) (purple dashed line), i.e. the geotherm dominated by a steep T gradient at 80 km depth, which occurs at the transition from partial to full coupling, as reported by Penniston-Dorland et al. (2015) and the average slab-top geotherm from Penniston-Dorland et al. (2015) (green dashed lines) are shown for comparison. Metamorphic facies modified from Peacock (1993) and Liou et al. (2004). Serpentine breakdown reactions from Padrón-Navarta et al. (2010) (PN10) and Hermann et al. (2000) (H00). Mineral abbreviations from Whitney and Evans  (2010).


During burial and heating the different slab lithologies undergo dehydration reactions. The produced fluid escapes from the source rock and migrates upward, likely interacting with the surrounding units of different chemical and isotopic composition. The effect of an external fluid input on the δ18O value of growing minerals is strongly dependent on the isotopic composition of the infiltrating fluid (δ18Ofluid) and on the degree of fluid–rock interaction. The integrated FR ratio is defined here as the total mass of aqueous fluid that has passed through and interacted with the rock normalized to the mass of the rock. To explore different scenarios, three models are discussed involving different associations of fresh or altered oceanic basalts with terrigenous or carbonate sediments (Fig. 1b): (1) the interaction between the fluid released from the metamorphosed mafic crust and the overlying metasediment is negligible and the two rocks evolve independently (no-interaction case, NI); (2) part of the fluid derived from the metabasalt (50 % when not specified differently) equilibrates with the metasediment, while the other part leaves the system (partial-interaction case, PI); and (3) all the fluid released by the metamorphosed mafic crust equilibrates with the metasediment before escaping the system (high-interaction case, HI). The fluid released by the entire system is a mixture of fluids derived from the progressive dehydration of the metabasalt and metasedimentary layers. In the NI case, both rock types behave like closed systems and the fluid is liberated from the metabasalt and from the metasediment separately. In the case of infiltration of fluid derived from the metabasalt in the metasediment, the amount of fluid released by this latter includes the fluid produced by dehydration reactions plus the excess fluid that enters the metasediment and cannot be incorporated into stable hydrous minerals.

The thickness and the degree of serpentinization of the lithospheric mantle subducting beneath the oceanic crust can be highly variable. The most important dehydration reactions in partly or fully serpentinized mantle are related to antigorite breakdown, which can release up to 12 wt % of water, playing an important role for water flows in subduction zones. Deserpentinization is assumed to result in two main subsequent fluid peaks (Padrón-Navarta et al., 2013; Scambelluri et al., 2004) related to the reactions


The effect on the δ18O of the metabasalts and metasediments of an external fluid influx, i.e. caused by dehydration of the underlying serpentinites, was investigated by defining an amount of fluid with a specific δ18O value that infiltrates the basaltic layer (Rock 1) at two steps of the model (480 and 660 C along the chosen PT path) (Figs. 1c, 2).

2.2 Model strategy

The strategy behind PTLoop consists of forward modelling the evolution of the mineral assemblage and the oxygen isotope composition of a rock column composed of two lithologies (Fig. 1a) of assigned thickness and starting bulk chemical and oxygen isotope compositions along a defined PT path using a stepwise procedure (Fig. 2). At each PT step, (1) the equilibrium mineral assemblage, oxygen isotope composition of stable phases (δ18O ‰ vs. Vienna Standard Mean Ocean Water, VSMOW), mass (in kg) and isotopic composition of the excess fluid for the metabasalt are calculated; (2) any fraction of excess fluid deriving from dehydration reactions in the metabasalt can be transferred to the metasediment or directly escapes the system; (3) the equilibrium mineral assemblage, δ18O value of stable phases, and amount and δ18O of the excess fluid for the metasediment are evaluated by accounting for the changes caused by the fluid input from the metabasalt; (4) the mass and δ18O of the total fluid leaving the system are calculated. Furthermore, at each step a chosen amount of external fluid with a given δ18O can be input into the metabasalt and its contribution is accounted for in the subsequent steps. This model is based on the assumption of thermodynamic equilibrium applied to a partially reactive system, whereby phases are assumed to reach chemical and isotopic equilibrium at all steps within the reactive part of the system, i.e. removing from the reactive bulk the phases that are fractionated (Lanari and Engi, 2017). Such petrological models can account for element sequestration during prograde metamorphism. Mineral fractionation in relics and fluid input or loss are two processes that are allowed to modify the reactive bulk composition. No mineral resorption is permitted. Any fluid liberated during dehydration (excess fluid) does not interact further and leaves the rock. This process is termed Rayleigh volatilization (Rumble, 1982; Valley, 1986). In natural rocks, it often occurs combined with the opposite endmember, the batch volatilization,where the produced fluid stays within the system as the mineral reaction proceeds and remains in isotopic equilibrium with the rock until the reaction is completed. In most natural cases involving oxygen isotopes, the difference between the results calculated using the two processes is negligible (Baumgartner and Valley, 2001). The released fluid is regarded as pure H2O. Any other effect related to solute transport by the fluid is ignored in the calculation; the potential effect of a CO2 component in the fluid is discussed in Sect. 3.2.

2.3 Governing equations

Equilibrium assemblage calculation for a given bulk rock composition at any P and T is performed with the program Theriak (de Capitani and Brown, 1987; de Capitani and Petrakakis, 2010) and is based on Gibbs energy minimization. A complete description of the minimization algorithm is given by de Capitani and Brown (1987). The reacting bulk composition may evolve in the course of the metamorphic history of a rock because of mineral fractionation, fluid loss or input of external fluids. The effective bulk composition is recalculated by PTLoop at each subsequent stage following the strategy of Lanari et al. (2017).

As for phase assemblage determination, equilibrium is a common assumption of stable isotope transport (Baumgartner and Rumble, 1988; Baumgartner and Valley, 2001; Bowman et al., 1994; Gerdes et al., 1995a, b). Thus, a molar equilibrium constant (K) can be defined to describe the thermodynamic stable isotope equilibrium between two substances i and j (Sharp, 2017):

(3) K = 18 O i / 16 O i 18 O j / 16 O j .

The fractionation factor (α) can be related to the equilibrium constant K as

(4) α = K 1 / n ,

where n is the number of exchanged atoms, normally 1 for simplicity. In isotope geochemistry, the isotopic composition is commonly expressed in terms of δ values:

(5) δ i = R i R St - 1 10 3 ( ) ,

where Ri and RSt are the isotope ratio measurements for the compound i and the defined isotope ratio of a standard sample respectively. For differences in δ values or for δ values of less than ∼10 ‰,

(6) 1000 ln α i - j = δ i - δ j

is a valid approximation that is used in most cases (Hoefs, 2015; Sharp, 2017). For oxygen isotope fractionation, the equation that can reproduce most of the available calibrations describing the stable isotope fractionation function between two phases is a second-order polynomial of 103T. Hence the stable isotope fractionation between two phases i (with k endmembers) and j (pure) as a function of T is described by Eq. (7):

(7) δ 18 O i - δ 18 O j = 1 k A k , j 10 6 T 2 + B k , j 10 3 T + C k , j X k , i N k , i N i .

The conservation of the bulk δ18O in the system is described by Eq. (8):

(8) δ 18 O sys N sys = k = 1 p M k N k δ 18 O k ,

where δ18Oi, δ18Oj and δ18Osys are the isotopic compositions of phase i, phase j and the system (bulk δ18O) respectively; Ak,j, Bk,j and Ck,j are the fractionation parameters for endmember k of mineral i vs. phase j; Xk,i is the fraction of endmember k in the phase i; Nk,i, Ni and Nsys are the total number of moles of oxygen in endmember k, in mineral i and in the system respectively; p is the number of phases; Mk is the number of moles of phase k; Nk is the its number of oxygen; and δ18Ok is its oxygen isotope composition. Given a stable mineral assemblage at any PT condition, the oxygen isotope partitioning among the stable phases is calculated by solving the linear system described by the sets of p−1 equatiomns of type (7) and Eq. (8); Kohn, 1993; Vho et al., 2020). In closed systems, the first term in Eq. (8) is constant. Open-system behaviour can either modify the δ18Osys or the number of moles of the phases (Nsys). The parameters A, B and C between phases were taken from the internally consistent database for oxygen isotope fractionation DBOxygen version 2.0.3 (Vho et al., 2020).

2.4 Starting assumptions

In order to represent the variability in the basaltic portion of the oceanic crust two different bulk compositions were used (Table 1): (1) a representative average MORB basalt (Gale et al., 2013) that has been hydrated but without any other addition or removal of element (metabasalt(h)) and (2) a basalt that underwent extensive sea-floor alteration during hydration (Baxter and Caddick, 2013 after Staudigel et al., 1996) that will be referred to as metabasalt(a) in the following. Those compositions are in good agreement with other compilations reported in the literature (e.g. Sun and McDonough, 1989; Albarède, 2005; Staudigel, 2014; White and Klein, 2014). Oceanic sediments were modelled with two distinct bulk compositions (Table 1): (1) terrigenous sediment referred to as metasediment(t) (clay from the Mariana Trench; Hacker, 2008, after Plank and Langmuir, 1998) and (2) nanno ooze carbonate referred to as metasediment(c) (Plank, 2014). Nanno oozes are widespread carbonate sea-floor sediments (Plank, 2014), and they are close in composition to carbonate-rich sediments observed in HP terrains (e.g. Bebout et al., 2013; Kuhn et al., 2005). Thicknesses of 1000 m for the basaltic layer, of 175 m for the clay sediment and 75 m for the carbonate sediment were chosen in order to maintain proportions between oceanic crust and sediments comparable with the values reported in various compilations (e.g. Hacker, 2008; Plank, 2014). This results in a total thickness of the rock column 2 to 3 times smaller than the real thickness to encompass the assumption of homogeneous temperature over the whole column within ∼20C. To overcome the effects of possible temperature variations within the column, a discretization step size of ∼20C along the PT path was applied.

Table 1Bulk compositions used for the simulations.

1 Gale et al. (2013). 2 Baxter and Caddick (2013) after Staudigel et al. (1996). 3 Mariana clay from Plank and Langmuir (1998). 4 Nanno ooze from Plank (2014). 5 Used for the thermodynamic modelling. 6 Water content at saturation at 350 C and 1.3 GPa. The calculated initial water content is consistent with values from the literature at these conditions (e.g. Poli and Schmidt, 1998; Hacker et al., 2003). NA: not available.

Download Print Version | Download XLSX

The bulk compositions were simplified to the Na2O–CaO–K2O–MgO–FeO–Al2O3TiO2SiO2H2O system; C is present in the initial bulk composition of the metabasalt(a) and the metasediment(c) (Table 1). MnO was excluded because it overemphasizes the stability of garnet at low metamorphic conditions (T≤350C). The conditions of the garnet-in reaction in Mn-absent systems match the results obtained for garnet nucleation in natural rocks (e.g. Laurent et al., 2018). Thermodynamic modelling was performed using the internally consistent dataset of Holland and Powell (1998) and subsequent updates (tc55, distributed with Theriak-Domino 4 February 2017; see Supplement S1). The following activity models where used for the solid solutions: Holland and Powell (2003) for calcite–dolomite–magnesite; Holland and Powell (1998) for garnet, white mica and talc; Holland and Powell (1996) for omphacite; Holland et al. (1998) for chlorite; Diener et al. (2007) for amphibole. In the presented model garnet undergoes fractional crystallization both in Rock 1 and Rock 2 fractionating from the reactive bulk for the subsequent steps. The amount of initial H2O in each rock was set at saturation and is reported in Table 1. No pore fluid expulsion, diagenetic and low-grade (T < 350 C) devolatilization reactions are considered in this study (see above).

A starting bulk δ18O for the metabasalt(h) of 5.7 ‰ was chosen and represents the reference value for an unaltered MORB (e.g. Cartwright and Barnicoat, 1999; Eiler, 2001; Staudigel, 2014; White and Klein, 2014), while a starting bulk δ18O of 9.0 ‰ is used for metabasalt(a), representative of basaltic material that underwent sea-floor alteration at T≤400C (e.g. Alt et al., 1986; Cartwright and Barnicoat, 1999; Eiler, 2001; Gregory and Taylor Jr., 1981; Miller and Cartwright, 2000; Staudigel, 2014; White and Klein, 2014). The starting bulk of the terrigenous sediment of 15 ‰ represents the average for the δ18O of clastic sediments reported by Eiler (2001). The chosen δ18O starting bulk of the carbonate sediment is 25 ‰, which represents a conservative estimate of marine carbonate δ18O (typically 25 ‰–35 ‰; Eiler, 2001). It is ∼5 ‰ higher than the values for metasedimentary carbonates in the Italian Alps (e.g. Cook-Kollars et al., 2014) that are likely to have interacted with lower- δ18O fluids during subduction.

In order to define the contribution of an external fluid originating in the lithospheric mantle by serpentine breakdown, a layer of 150 m of pure serpentine containing 12 wt % bulk H2O was considered and the mass of water released at each reaction was calculated by mass balance, resulting in an input of 7800 kg of water at 480 C and of 25 350 kg at 660 C to satisfy the reactions in Eqs. (1) and (2) respectively. In order to fit the thicknesses chosen for the oceanic crust and the sedimentary layer (2 to 3 times thinner than an average lithospheric section), the 150 m of pure serpentine correspond to a conservative estimate of 3000 m of serpentinized peridotite with an average serpentine content of 5 % in volume. This is in agreement with the values used by Barnes and Straub (2010) and John et al. (2011) based on the estimate by Sharp and Barnes (2004). Serpentine oxygen isotope compositions reported in the literature are highly variable (Cartwright and Barnicoat, 1999, 2003; Früh-Green et al., 2001; Mével, 2003; Miller et al., 2001), typically ranging from 1 to 10 ‰. In mid-oceanic ridge environments, the distribution has a peak between 2 ‰ and 5 ‰ (Mével, 2003). A value of 2.5 ‰ was chosen from the lower-δ18O side of this peak. This results in the liberation of a fluid with a characteristic low-δ18O value, but one which is still feasible for natural serpentinites; this value is clearly distinct with respect to the overlying lithologies. The effect of infiltration into the metabasalts and metasediments of serpentinite-derived fluids will become smaller as the fluid δ18O becomes higher, approaching the equilibrium with the overlying lithologies. The δ18O value of the released fluid is ∼4.5 ‰ at T>550C (serpentine–water oxygen isotope fractionation factors from Vho et al., 2020). Further details on the modelling input data are given in Supplement S2.

3 Results

3.1 Stable mineral assemblage

The evolving stable mineral assemblages and bulk water contents of each lithology, without external fluid input, were calculated for each rock composition along the prograde PT path. Results are provided as mode-box diagrams in Fig. 3. The H2O field represents the volume fraction of excess water in each rock type. The fluid is progressively extracted becoming isolated from the reactive part of the system. Garnet is the only phase prevented from re-equilibrating in the model, thus fractionating from the reactive bulk composition. Below 450 C and 1.80 GPa, in the metabasalts glaucophane, actinolite and lawsonite comprise ∼80 vol. % of the paragenesis, with minor phengite (Siapfu=3.67–3.63, XMg=0.62–0.56 in metabasalt(h); Siapfu=3.68–3.65, XMg=0.67–0.58 in metabasalt(a)), omphacite (XNa=0.45–0.42, XMg=0.81–0.72), chlorite and titanite. Metabasalt(h) is richer in SiO2, FeO and MgO with respect to the metabasalt(a), and chlorite is stable up to 480 C. Metabasalt(a) contains ∼5 vol. % of Ca carbonate that remains stable over the entire PT path. For either composition, the volume of glaucophane, actinolite and lawsonite gradually decreases from 480 C and ∼1.90 GPa until complete consumption at 600–620 C and 2.30–2.36 GPa. Those represent the major hydrous phases contributing to the dehydration, while a secondary role is played by talc and zoisite at higher PT conditions (T≥580C, P≥2.24 GPa). Most of the water still retained in the rocks is stored in phengite, the abundance of which is primarily controlled by bulk K2O content, higher in the metabasalt(a), and that remains stable beyond the model conditions. Garnet production starts at ∼500C and ∼2.00 GPa (Xalm=0.60, Xgrs=0.35), and garnet grows continuously until constituting ∼20 vol. % of the metabasalt(a) and ∼35 % of the metabasalt(h).

Figure 3Mode-box diagrams showing the evolution of the mineral assemblages and fluid during subduction of the different rock types along the geotherm shown in Fig. 2. The initial H2O content is < 1 vol % in the metabasalts and in the metasediment(c) and ∼3 vol. % in the metasediment(t) (see Table 1 for details). Garnet fractionation is applied to all the lithologies. The volume fraction of garnet shown at each step represents the sum of the fractionated and newly grown garnet. The phase proportions refer to the NI case, where the H2O content is the excess H2O produced by each rock type evolving independently. The excess water is fractionated at each step and the volume fraction displayed represents the sum of the fractionated and the newly produced water. Mineral abbreviations from Whitney and Evans (2010).


In the metasediment(c) calcium carbonate, quartz, phengite (Siapfu=3.42, XMg=0.53) and jadeite (jadeite content Xjd=0.95, XNa=0.96, XMg=0.39) compose ∼80 vol. % of the solids at 350 C and 1.30 GPa. Glaucophane and lawsonite are stable up to 460 and 560 C respectively. Jadeite abundance increases from ∼10 vol. % at 440 C to ∼16 vol. % at 460 C; ankerite is stable in minor amounts (≤3 vol. %) at 440–560 C. Garnet (Xalm=0.57, Xgrs=0.41) is stable only at 540–580 C and 2.12–2.24 GPa, reaching ∼5 vol. %, and is then preserved because of the assumption of fractionation from the bulk in the model. The metasediment(t) at T<500C shows a paragenesis of phengite (Siapfu=3.68–3.55, XMg=0.53–0.44), glaucophane, lawsonite, quartz and omphacite (XNa=0.50–0.48, XMg=0.87–0.70), with minor titanite. At 520 C and 2.06 GPa, lawsonite is consumed and the amphibole proportion reduces from ∼30 vol. % to < 20 vol. %, while garnet (Xalm=0.63, Xgrs=0.34) is produced and reaches ∼10 vol. %. At these conditions, the omphacite content decreases and a clinopyroxene of more jadeitic composition (Xjd=0.72, XNa=0.76, XMg=0.45) becomes stable. These models are in line with the first-order mineralogical changes observed in subducted (and exhumed) crustal material. Thermodynamic calculations predict the coexistence of a calcic amphibole and a sodic amphibole in the metabasalts and of jadeite and omphacite in the metasediment(t). From an oxygen isotope partitioning perspective, the interpretation of the modelled coexistence of a sodic and a calcic amphibole either as two endmembers of a solid solution or as coexisting minerals is equivalent. The same applies to pyroxenes, for which the modelled coexistence of two pyroxenes can be interpreted as a continuous solid solution. Therefore, this does not affect the oxygen isotope partition model final results for the other phases and the bulk.

3.2 Production of aqueous fluid

At the initial conditions, all the lithologies are saturated in H2O (Table 1). Up to 500 C, lawsonite, actinolite and glaucophane are the main repositories of H2O in the metabasalts, followed by chlorite and minor phengite. A significant pulse of water is modelled at 500–520 C and 2.00–2.06 GPa for both the metabasalts (Fig. 3a, b). This pulse is caused by a decreasing abundance of lawsonite and amphibole and a breakdown of chlorite followed by growth of garnet and omphacite. This first dehydration stage releases ∼25 % of the total water loss from the metabasalt(a) (∼4.0 vol. % H2O liberated) and ∼45 % from the metabasalt(h) (∼6.5 vol. % H2O liberated). The second significant pulse in the metabasalts occurs at 620–640 C and 2.36–2.42 GPa, releasing ∼40 % of the total water loss from the metabasalt(a) and ∼15 % from the metabasalt(h). This pulse is caused by the final breakdown of lawsonite and of amphibole in the metabasalt(a) (Fig. 3b). In the metabasalt(h), glaucophane and actinolite breakdown takes place at 600 C and 2.30 GPa. Growth of talc only incorporates half of the water released by amphibole breakdown, causing an intermediate fluid pulse of minor magnitude at these PT conditions (Fig. 3a). If an initial undersaturated basaltic composition was considered, the amount of released fluid due to the breakdown of hydrous phases would be smaller.

The metasediment(c) is the rock type that dehydrates the least: the two main pulses of fluid production are at 480 C and 1.92 GPa (∼0.4 vol. % H2O liberated) and from 540 C and 2.12 GPa to 560 C and 2.18 GPa (∼1.7 vol. % H2O liberated), caused by a breakdown of glaucophane and lawsonite respectively (Fig. 3c). The water produced from these two dehydration stages represents < 0.02 wt % of the total water released by the system composed of metabasalt(h) and metasediment(c). In the metasediment(t), the main fluid pulse occurs at 520 C and 2.06 GPa (∼3.0 vol. % of H2O liberated), caused by the breakdown of lawsonite and a decrease in glaucophane abundance (Fig. 3d). The water produced from this dehydration stage represents ∼0.07 wt % of the total water released by the system composed of metabasalt(a) and metasediment(t).

As specified in Sect. 2.2, the aqueous fluid used for the calculation considers only the water component (aH2O=1). The release of CO2 occurs during progressive metamorphism of metabasalt(a) and metasediment(c). CO2 is absent or present in negligible amounts (< 1 mol % of the total fluid) up to 440 C and 1.76 GPa. At higher temperatures, the X(CO2) content increases significantly up to ∼10 mol % in the fluid released by metabasalt(a) and ∼30 mol % in the fluid released by metasediment(c) at 700 C and 2.60 GPa. However, the total amount of fluid produced at these conditions is negligible (< 0.01 vol. %). The oxygen isotope fractionation database used for this study (Vho et al., 2020) does not include fractionation data for CO2. Limited calibrations are available for oxygen isotope fractionation between H2O and CO2 (Friedman and O'Neil, 1977; O'Neil and Adami, 1969; Zheng, 1994). The fractionation values diverge significantly, for example between 4.9 ‰ and 8.3 ‰ at 450 C and between 1.5 ‰ and 4.4 ‰ at 700 C. At 520 C, the fluid released by metabasalt(a) contains 2 mol % of CO2 and 6 mol %–7 mol % at 620–640 C. The consideration of CO2 would produce a negligible shift on the fluid δ18O (0.1 ‰–0.2 ‰ at 520 C and fractionation used). The amount of CO2 in metasediment(c) derived-fluid is 2 mol % at 480 C and 5 mol % at 540 C. The consideration of CO2 would also produce negligible to minor shifts in the fluid δ18O (0.1 ‰–0.3 ‰ at 480 C and 0.2 ‰–1.1 ‰ at 620 C, depending on the fractionation).

3.3 Oxygen isotope composition

The largest initial bulk δ18O difference occurs between metabasalt(h) and metasediment(c) (14.3 ‰, the relatively water-poor system), while the smallest initial bulk δ18O difference is observed between metabasalt(a) and metasediment(t) (6.0 ‰, the relatively water-rich system). In the following, the results are presented in detail for a selection of two endmember scenarios (Figs. 4, 5): (1) metasediment(c) associated with metabasalt(h) and (2) metasediment(t) associated with metabasalt(a) when not specified differently. Other scenarios (i.e. metabasalt(h) associated with metasediment(t) and metabasalt(a) associated with metasediment(c)) give intermediate results in terms of oxygen isotope composition variations as a consequence of fluid–rock interaction. Further details and the results for the intermediate scenarios are given in Supplement S3.

3.3.1 Bulk oxygen isotope compositions

For rocks that undergo only dehydration reactions, the starting bulk δ18O evolves as a consequence of garnet and fluid fractionation (Fig. 4). The bulk δ18O shift related to water fractionation is within 0.2 ‰, while the shift due to garnet fractionation is within 0.5 ‰ (Fig. 4c, d). Since water has typically heavier δ18O with respect to the bulk and the garnet a lighter one, the two effects produce opposite trends. The combination of both effects results in a shift in the bulk δ18O in the considered lithologies restricted to < 0.3 ‰. This in turn leads to negligible (< 0.2 ‰) variations in the δ18O values of the stable phases.

Figure 4Calculated bulk and mineral δ18O values along the geotherm shown in Fig. 2. Bulk δ18O: black solid line. Hydrous mineral δ18O: coloured dashed lines. Anhydrous mineral δ18O: coloured solid lines. Released H2O: thick blue dashed lines. (a, b) Modelled mineral, bulk and released fluid δ18O values for the metabasalt(h) and the metabasalt(a) considering garnet fractionation and excess fluid loss and in the absence of external fluid input. (c, d) Quantification of the effects of garnet fractionation and fluid loss on the bulk δ18O of the metabasalts. (e, f) Modelled mineral, bulk and released fluid δ18O values for the metasediment(c) and the metasediment(t) considering garnet fractionation and excess fluid loss and in the absence of external fluid input. Scale on the top x axis indicates P (GPa). Mineral abbreviations from Whitney and Evans (2010).


In the metasediments, the progressive interaction with the fluid from the metabasalts causes a decrease in the bulk δ18O (Fig. 5) that is controlled by the amount and δ18O of the incoming fluid. A significant decrease starts at 480 C, where the amount of water released by the metabasalts increases of about 1 order of magnitude from < 0.05 vol. % to ∼0.3 vol. % due to partial consumption of amphibole and lawsonite. The maximum shift in bulk δ18O was calculated for the metasediment(c) interacting with the metabasalt(h)-derived fluid at 12.9 ‰ for the HI case (integrated FR ratio = 0.75 kg kg−1), while it is 8.7 ‰ for the PI case (integrated FR ratio = 0.38). The shift in the bulk δ18O of the metasediment(t) interacting with the metabasalt(a)-derived fluid is 1.5 ‰ for the HI case (integrated FR ratio = 0.35) and 2.7 ‰ for the PI case (integrated FR ratio = 0.18).

Figure 5Calculated bulk and mineral δ18O values in the metasediments along the geotherm shown in Fig. 2 in the case of interaction with the metabasalt-derived fluid. Only bulk, released H2O and representative mineral δ18O values are shown for clarity. Garnet fractionation and excess fluid loss are considered. (a) Metasediment(c): evolution of the δ18O values in the case of NI (continuous lines) and PI (dashed lines). (b) Metasediment(t): evolution of the δ18O values in the case of NI (continuous lines) and PI (dashed lines). (c) Metasediment(c): evolution of the δ18O values in the case of NI (continuous lines) and HI (dashed lines). (d) Metasediment(t): evolution of the δ18O values in the case of NI (continuous lines) and HI (dashed lines). Scale on the top x axis indicates P (GPa). Mineral abbreviations are from Whitney and Evans (2010).


3.3.2 Oxygen isotope composition of mineral phases

Since at infinite temperature the fractionation between any two phases approaches 0 ‰, a general trend in the reduction of oxygen isotope fractionation between the stable phases with increasing metamorphic grade is observed in all lithologies and is a result of the temperature increase (Fig. 4a, b, e, f). As a consequence, mineral phases typically heavier than the bulk (i.e. quartz and carbonates) become isotopically lighter with increasing metamorphic conditions, and the mineral phases typically lighter than the bulk (i.e. rutile, garnet and titanite) become isotopically heavier. Such variations are limited (i.e. within 1.0 ‰) for most of the phases, with the exception of quartz, calcite and rutile that may vary up to 3.0 ‰ in response to temperature variation only in the considered range.

In the case of ingress of the low-δ18O fluid from the metabasalts in the metasedimentary rocks (PI and HI cases), the mineral δ18O values decrease progressively with respect to the NI case following the trend of the bulk δ18O (Fig. 5). For instance, in the case of NI the δ18O of quartz in the metasediment(t) decreases from 19.4 ‰ to 17.3 ‰ (2.1 ‰) and that of quartz in the metasediment(c) from 28.0 ‰ to 26.5 ‰ (1.5 ‰) over the total temperature range modelled. In the PI and HI cases, the δ18O shift in quartz is respectively 3.8 ‰ and 5.0 ‰ in the metasediment(t) and 10.0 ‰ and 13.9 ‰ in the metasediment(c). Hence, the final quartz δ18O values (at 700 C, 2.60 GPa) for the HI case in the metasediment(t) and in the metasediment(c) are respectively 3.0 ‰ and 13.0 ‰ isotopically lower than the expected values in the case of NI (Fig. 5). The maximum shift in δ18O (i.e. between NI and HI cases) for the other stable phases is within those values. In the metasediment(t), the δ18O values of phengite, glaucophane, jadeite, rutile and garnet decrease by 3.0 ‰ and those of omphacite by 2.4 ‰ from the NI to the HI case. Lawsonite and titanite are not stable after the first significant input of metabasalt-derived fluid (500 C, 2.00 GPa), and their δ18O values decrease by 1.1 ‰ and 0.3 ‰ respectively from the NI to the HI case. In the metasediment(c), the δ18O values of dolomite, jadeite, phengite, rutile and aragonite decrease by a maximum of 13.1 ‰, of lawsonite and ankerite by a maximum of 10.3 ‰ and 9.5 ‰ respectively. Garnet crystallizes only between 540 and 580 C and its δ18O value in the HI case is 5.9 ‰ lower than the one in the NI case.

3.3.3 Oxygen isotope composition of the fractionated fluids

The δ18O of the fractionated water from each rock type at each step is in isotopic equilibrium with the stable mineral assemblage at the given conditions. In the temperature range where most of the fluid is released (i.e. T>480C), the δ18O of the fluid is 7.0±0.5 ‰ for the metabasalt(h) and 10.0±0.5 ‰ for the metabasalt(a). At T>480C, the water released in the NI case by the metasediment(c) has a δ18O= 24.2–25.6 ‰ and that released from the metasediment(t) has a value of 15.4 ‰–16.4 ‰ (Fig. 6). The δ18O of the fluid leaving the system, e.g. infiltrating an upper layer or the mantle wedge, results from the mixing of the water released from the metabasalts and the overlying metasediments, and it derives from the balance between the amount of fluid released by each rock type and its δ18O value (Fig. 6a, b).

Figure 6Double plot diagrams showing the oxygen isotope composition of the released fluids (left axis, coloured lines) and the amount (in wt % of the total rock) of the total fluid released by the systems (right axis, grey field) for each interaction case in the absence of serpentinite-derived fluid input. (a) Modelled fluid δ18O values and amount for the system metabasalt(h)+ metasediment(c). (b) Modelled fluid δ18O values and amount for the system metabasalt(a)+ metasediment(t). Dashed lines show the δ18O values of the fluids released by each rock type and solid lines the δ18O values of the final fluids released by each system. In the case of HI, all the metabasalt-derived fluid infiltrates the metasediment. Hence the final fluid released overlaps with the fluid expelled by the metasediment and only one line is represented (light blue, marked as HI). Because the δ18O values of the fluids released by the metabasalts are not affected by the degree of interaction, all three cases are represented by one line (marked as NI).


In the NI case, the δ18O of water leaving the system is up to 1 ‰ higher than the composition of the fluid released by the metabasalt because of the minor input from the metasediment at around 500–550 C. The only exceptions are for the interaction between the metabasalt(a) and the metasediment(t) at T of ∼450 and ∼700C. At these conditions, the δ18O values of the final fluid are up to 5 ‰ higher than the metabasalt(h)-derived fluid (Fig. 6b). This increase is caused by a predominance of the metasediment-derived fluid at those conditions; however, the amount of high-δ18O fluid represents < 0.1 wt % of the rock column and ∼1 wt % of the total released fluid (Fig. 6b).

In the case of interaction between the metasediments and the metabasalt-derived fluid, part of or all this fluid reacts with the overlying metasediment before leaving the system. The final δ18O of the fluid is controlled by the integrated FR ratio in the metasediments and their buffering capacity. In the HI case, the δ18O of the released fluid has a dominant sedimentary signature at T < 500–520 C, before the first fluid pulse from the metabasalt (14.5 ‰–15.5 ‰ for the metasediment(t) and 23.0 ‰–24.0 ‰ for the metasediment(c), Fig. 6). The first metabasalt-derived fluid pulse (500–520 C, see above) causes a drop in the bulk δ18O values of the total released fluids of 0.7 ‰ for the metabasalt(a)–metasediment(t) association and of 6.3 ‰ for the metabasalt(h)–metasediment(c) association. The second metabasalt-derived fluid pulse (620–640 C, see above) causes a second decrease in the δ18O values of the total released fluid equal to 1.0 ‰ for both the lithological associations.

3.3.4 Input of serpentinite-derived fluid

Ultramafic rocks tend to undergo episodic dehydration (see above). In the following, the effects caused by the input of a serpentinite-derived fluid at the base of the rock column in the case of HI are described (Fig. 7). The input of the amount of water with a δ18O of 4.5 ‰ (see above), corresponding to the dehydration of 150 m of pure serpentine at 480 C (2.0 wt % H2O) and 660 C (6.5 wt % H2O), has a limited impact on the δ18O of the system. It produces a decrease of < 0.2 ‰ in the final bulk δ18O of the metabasalts, up to ∼1.0 ‰ in the metasediment(c) (Fig. 7a) and < 0.5 ‰ in the metasediment(t) with respect to the HI case with no serpentinite-derived fluid (Fig. 7b). The largest decrease occurs at 660 C, where the second pulse of serpentinite-derived fluid enters the system, resulting in ca. 0.1 ‰, 0.3 ‰ and 1.0 ‰ for the metabasalts, metasediment(t) and metasediment(c)δ18O bulk values respectively (Fig. 7). Even by increasing the thickness of the serpentinite by a factor of 2, the variations in bulk δ18O values with respect to the HI case with no serpentinite-derived fluid are < 1.0 ‰ for any rock type with the exception of the metasediment(c), for which the bulk δ18O decreases by 2.0 ‰ with respect to the HI case with no serpentinite-derived fluid (Supplement S3). The effect of the serpentinite-derived fluid on the metabasalts remains the same for any interaction case (NI, PI and HI), while the variations in the δ18O of metasedimentary rocks decreases to zero with decreasing infiltration of external fluid (Supplement S3).

4 Discussion

4.1 Effect of stable assemblage evolution and phase fractionation on the δ18O

The changes in mineral assemblage, modes and compositions along a prograde PT path control (1) the oxygen isotope partitioning between the stable phases and (2) the amount and δ18O of the fluid released by the system. At the same time, oxygen isotope fractionation between the stable phases is controlled by temperature. Thus, the effects of evolving paragenesis and increasing temperature are systematically overlapping in nature. In the case of a closed system, the bulk concentrations of 18O and 16O remain constant and a change in one phase is compensated for exactly by adjustments in other phases (Baumgartner and Valley, 2001; Kohn, 1993). In this situation, major changes in mineral assemblage do not play a significant role in shifting the δ18O of stable phases: this is demonstrated by the limited (< 0.5 ‰) shift in δ18O values of quartz, garnet, phengite, omphacite and rutile in the metabasalt(h) after (1) the breakdown of amphibole and lawsonite and (2) the crystallization of talc and kyanite over a narrow temperature range between 500 and 580 C (Fig. 4a, b). Similar effects can also be anticipated for rocks with different chemical composition that undergo major changes in the mineral assemblage (see Supplement S4).

In a closed system evolving at equilibrium, the initial chemical bulk composition and bulk δ18O do not change along the PT evolution. However, in metamorphic rocks, compositional zoning and metamorphic overgrowths are often preserved in refractory minerals (Lanari and Engi, 2017) indicating that parts of the minerals have become isolated from the reactive volume of the rock. This scenario is commonly referred to as “partially re-equilibrated open systems”, because the chemical and the isotopic compositions vary as a consequence of the fractionation of solid and fluid phases (i.e. garnet fractionation and excess fluid removal) even in the absence of external fluid input. Phase fractionation is expected to affect the bulk δ18O as function of both the amount of fractionated or expelled phases and their isotopic composition. Fractionation of a phase lighter than the bulk in δ18O leads to an increase in the reactive bulk δ18O value, while fractionation of an isotopic heavier phase leads to a decrease in the reactive bulk δ18O value.

The most common example of fractionating metamorphic mineral is garnet, which systematically records compositional zonation at low- to medium-grade (Evans, 2004; Giuntoli et al., 2018; Konrad-Schmolke et al., 2008; Spear, 1988; Tracy, 1982). Therefore, garnet fractionation was incorporated into the model in order to better approximate the behaviour of natural systems. Note that this effect is reduced at higher grade where intra-crystalline diffusion becomes efficient to partially re-equilibrate garnet (Caddick et al., 2010; Lanari and Duesterhoeft, 2019). As already documented by Konrad-Schmolke et al. (2008), garnet fractionation controls the extent of the garnet stability field. Garnet crystallization is not systematically expected to occur near the peak conditions, if the matrix was strongly depleted due to garnet fractionation and the volume of garnet remains constant (i.e. for the metabasalt(a), Fig. 3b). While garnet fractionation is recognized to significantly affect isopleth thermobarometry and volume fractions (Lanari and Engi, 2017), its effect on oxygen isotope bulk composition and partitioning is negligible (< 0.5 ‰) in all the studied lithologies. In the model, the garnet fraction varies from ∼5 vol. % in the metasediment(c) to ∼35 vol % in the metabasalt(h) (Fig. 3), and its δ18O is 0.8 ‰ to 1.7 ‰ lower than the bulk (Fig. 4a, b, e, f).

Figure 7Double plot diagrams showing the effect of the fluid deriving from a layer of 150 m of pure serpentine (see text for details) on the δ18O of the rock types and of the of the total released H2O (left axis) and on the amount and distribution of the H2O released by the systems (right axis). All the values are calculated assuming HI between the metabasalts and the metasediments. Black dashed lines represent the bulk δ18O of the different lithologies and the blue lines the δ18O of the final fluid released by the systems. The final H2O released by each system is represented with a red line. The amount and distribution of the final fluids in the case of no serpentinite-derived fluid input (grey fields) are shown for comparison. (a) Modelled δ18O values and fluid amount for the system metabasalt(h)+ metasediment(c). (b) Modelled δ18O values and fluid amount for the system metabasalt(a)+ metasediment(t).


Beside garnet fractionation, dehydration due to hydrous mineral breakdown and expulsion of excess water may contribute to changing the starting chemical and isotopic bulk compositions. Baumgartner and Valley (2001) postulated that the liberation of metamorphic fluids might have a profound effect on the stable isotope composition of the residual rock. In the present study, the maximum fluid loss is from the metabasalts that release ∼15 vol. % (∼5 wt %) of H2O with δ18O values 0.9 ‰ to 1.5 ‰ higher than the bulk rock (Fig. 4a, b) at T≥500C. This significant fluid flux produces a decrease in the bulk δ18O of less than 0.2 ‰ (Fig. 4c, d). Even if more extensive dehydration occurs, the effect on the bulk δ18O value will be typically lower than 1.0 ‰. No significant differences in the effect of stable assemblage evolution and phase fractionation are observed between the four lithologies. Therefore, the bulk δ18O of a rock that experienced a succession of dehydration reactions, without rehydration by external fluids or major compositional changes through decarbonation or mineral dissolution, is likely to be representative of its protolith composition. In this regard, integrated thermodynamics and oxygen isotope modelling represent a key tool for quantifying the potential effects of different processes and for assessing closed- or open-system behaviours.

4.2 Mineral δ18O zoning as indicator of open-system behaviour

In the last decades, the significant advances of oxygen stable isotope analyses by SIMS (secondary ion mass spectrometer) have allowed zoned metamorphic minerals to be analysed in situ with a precision down to 0.2 ‰–0.3 ‰ (2σ) (e.g. Ferry et al., 2014; Martin et al., 2014; Page et al., 2010). The magnitude of the intra-crystalline δ18O variation in key metamorphic minerals has been widely used to establish whether a metasomatic stage is related to an internal fluid deriving from the breakdown of hydrous phases or if it reflects equilibration with an external fluid of different isotopic composition (e.g. Putlitz et al., 2000; Errico et al., 2013; Page et al., 2013; Russell et al., 2013; Martin et al., 2014; Rubatto and Angiboust, 2015; Engi et al., 2018). Understanding the scale of fluid migration at depth and the magnitude of the interaction between fluids and minerals is of special interest and can be enhanced by modelling of such fluid flow and isotopic exchange (Baumgartner and Valley, 2001). The definition of different interaction cases (NI, PI, HI) is useful to represent various degrees of isotopic exchange between the fluid and the rock. If the flow is channelled, all the interaction cases can possibly occur in close proximity and the modelled scenarios can be linked to the evolution of different domains around the vein (Fig. 8). The flow of a pervasive fluid leads to the homogenization of the chemical potential of all components, including stable isotopes (Baumgartner and Valley, 2001), and it is represented by the HI case, as long as integrated FR ratios are high. In contrast, the flow of a channelled fluid results in local chemical heterogeneities, allowing some portions of the rock and of the fluid to remain unaffected (NI case) and some others to be only partially affected (PI case).

Figure 8Schematic section of a channelled fluid flow where different degrees of exchange between the fluid and the rock may occur in spatial proximity. From the host rock perspective, the NI case describes the distal portion of the rock walls where no fluid infiltrates, the PI case the intermediate portion where a limited amount of external fluid is available and the HI case the pervasively infiltrated rock proximal to the vein. From the fluid perspective, the NI case describes the fluid flow in the centre of the channel, for which the exchange with the rock walls is negligible; the PI case the situation in which part of the fluid does not react with the wall rock and part equilibrates with it, and the HI case the situation in which no fluid flow occurs without equilibration with the host rock.


The first step for a meaningful interpretation of an observed intra-grain variation in δ18O value is the quantification of the possible effects of changes in T and mineral assemblage. Such effects are characteristic of each phase (Fig. 4a, b, e, f). Quartz, calcite and rutile are the minerals most sensitive to temperature changes. Their composition is expected to vary up to 1 ‰ per 100 C, and they are stable over a wide range of temperature. For such phases, care is required in interpreting significant intra-grain δ18O variations (i.e. up to 3.0 ‰) since it does not necessarily reflect interaction with an external fluid having a different isotopic composition. However, quartz and calcite are not necessarily robust minerals that can preserve chemical and isotopic zoning during metamorphism. On the other hand, the variation in δ18O value over 150 C in a mineral such as garnet that commonly retains growth zoning is typically within 0.5 ‰ when no phase fractionation is involved and still less than 1.0 ‰ when considering the effect of previous garnet and/or excess fluid fractionation (Fig. 4a, b, e, f). Thus, any larger variation has to be linked to a significant change in bulk δ18O. Similar behaviour is observed for other key metamorphic minerals such as phengite, amphibole and clinopyroxene. These minerals have been widely used in metamorphic petrology as thermometers and geobarometers (Dubacq et al., 2010; Ferry and Spear, 1978; Parra et al., 2002) and are expected to be robust targets to link the fluid evolution along the PT path, especially when mineral relics are preserved. Due to its large capacity to preserve growth chemistry, garnet has been a primary target for microscale measurement of oxygen isotopes. Protocols and reference materials for SIMS measurements for a range of garnet compositions are well established (e.g. Martin et al., 2014; Page et al., 2010; Vielzeuf et al., 2005b) and its retentivity to high T resetting by diffusion has been investigated (Higashino et al., 2018; Vielzeuf et al., 2005a). In HP rocks, various degrees of intra-grain variations in garnet δ18O associated with external fluid infiltration have been reported in the literature. Where field constraints on the fluid source (i.e. oxygen isotope measurements on the feasible source lithologies) are available, the intra-grain variation can assist in the calculation of FR ratios. Martin et al. (2014) describe a shift of 2.5 ‰ associated with infiltration of serpentinite- or altered gabbro-derived fluids in metasediments from the continental basement in Corsica. Rubatto and Angiboust (2015) observed a shift of 3.5 ‰–4.0 ‰ in garnet from an eclogite from Monviso that they attributed to sediment-derived fluid infiltration. Vielzeuf et al. (2005b) measured a decrease of 2.5 ‰–3.0 ‰ between garnet core and rim in the Dora Maira gneiss probably caused by the interaction with a fluid derived from the whiteschists (Gauthiez-Putallaz et al., 2016). Studies that used in situ measurement of micas are limited (Bulle et al., 2019; Siron et al., 2017), and thus the potential of mica to trace fluid–rock interaction is still underexplored. The matrix complexity of pyroxenes and amphiboles remains a challenge for SIMS measurements.

4.3 Interaction with fluids from an ultramafic source

The effect of pervasive fluid flow deriving from the breakdown of serpentine (δ18O= 4.5 ‰) in a serpentine layer of different thicknesses (150, 300, 600 m, see above) on the bulk δ18O of the two metabasalts is negligible (< 0.5 ‰, Fig. 7, Supplement S3). This is mainly due to (1) the minor difference in δ18O between the serpentinite-derived fluid and the metabasalts (1.2 ‰ for the metabasalt(h) and 4.5 ‰ for the metabasalt(a)) and (2) the very low integrated FR ratio (0.01, 0.02 and 0.04 for the three cases). In this case, the integrated FR ratio seems to play a bigger role, since the limited change in δ18O is similar for both the metabasalts even if the initial Δ18O between fluid and rock is larger for the metabasalt(a). With the same total volume of fluid and rock, a channelled fluid flow would imply larger volumes of fluid interacting with smaller volumes of rock (higher FR ratios) and would thus be expected to drive larger variations in isotopic composition. For instance, by increasing the FR ratio by a factor of 10 (from 0.01 to 0.1) the bulk δ18O decreases by 0.6 ‰ (metabasalt(h)) and 1.1 ‰ (metabasalt(a)) upon infiltration of the serpentine-derived fluids.

In contrast to the metabasaltic layers, the overlying metasediments have (1) larger compositional differences with the serpentinite-derived fluid and (2) a smaller mass. The effect of the serpentinite-derived fluid input at the base of the column on the metasediment bulk δ18O compositions can be up to 10 times larger than the effect that the same amount of fluid has on basaltic compositions. This is the case even when the serpentinite-derived fluid completely equilibrates with the metabasalt before interacting with the overlying metasediment. In our model, the fluid infiltrating the metasediment is always a mixture of metabasalt-derived and serpentinite-derived fluid. Considering instead an input of serpentinite-derived fluid directly in the metasediment and applying an FR ratio of 0.1, the final bulk δ18O of the metasediment(c) decreases by 4.6 ‰ and that of the metasediment(t) by 2.6 ‰. These values are significant and are comparable with variations observed in natural samples. For instance, Martin et al. (2014) describe a shift in δ18O of 2.5 ‰ among different generations of HP garnet in a sample from the Corsica continental basement (garnet mantle δ18O=7.2±0.4 ‰, garnet rim δ18O=4.7±0.5 ‰). The authors associate this shift with an infiltration of serpentinite-derived fluids and, to a lesser extent, altered gabbro-derived fluid. Williams (2019) describe an extreme δ18O shift of 15 ‰ between garnet core and rim in a metasediment from the Lago di Cignana Unit. Such an oxygen isotope composition variation has been related to a strongly channelled fluid influx originating from the dehydration of serpentinites. These results demonstrate that metasediments can be a good target to detect fluids from an ultramafic source migrating upward through the subducting slab or along the subduction interface, even though the two lithologies may not be in direct contact in the field.

It is important to note that the relatively water-poor system composed of the metabasalt(h) and the metasediment(c) is more sensitive to external fluid infiltration and thus affected by the highest changes in δ18O, according to observations in natural systems (e.g. Page et al., 2019).

4.4 Effect of the subduction geotherm

As discussed in detail by previous studies (Baxter and Caddick, 2013; Hacker, 2008; Hernández-Uribe and Palin, 2019; Syracuse et al., 2010), the subduction geotherm has an important effect on hydrous phase stability and PT conditions of fluid release into the mantle wedge. Along the average D80 geotherm by Syracuse et al. (2010) (Fig. 2), the top of the slab crust releases ∼95 % of the water at ∼80 km, during the transition from partial to full coupling. Along the cold geotherm by Penniston-Dorland et al. (2015), the first significant fluid pulse (20 %–40 % of water released) occurs as a consequence of the breakdown of glaucophane and actinolite at greater depths than predicted by our model (∼500C and 2.7–2.8 GPa, 90–100 km depth), and the remaining water is released at depth > 100 km. Both those models imply a relatively dry mantle forearc region, contradicting what has been described by Bostock et al. (2002). By contrast, along the warm geotherm by Penniston-Dorland et al. (2015), the breakdown of chlorite, epidote and actinolite releases 40 %–50 % of the water at 460–470 C and ∼0.6 GPa (∼20 km depth); after this stage, ∼95 % of the water is released at depth < 70 km. This implies that little to no water is available at subarc depths. The average geotherm by Penniston-Dorland et al. (2015) is hotter than the one by Gerya et al. (2002) used in our model for T<650C and P<2.5 GPa (Fig. 2). Along this PT gradient, chlorite is the main water carrier in low-grade conditions and the first two main dehydration pulses occur at depths of 30–40 km (significant decrease in chlorite, 20 %–25 % of water released) and of ∼50 km (complete breakdown of chlorite, 20 %–25 % of water released). Along this geotherm, the most water is released at shallower depths than in our model. Differences can be investigated in detail by modelling each case with PTLoop. Nevertheless, the effects of fluid–rock interaction on the bulk and mineral δ18O compositions follow the general trends described above. Different parageneses are expected to form during hydration, but the shift in bulk δ18O remains constrained by the FR ratio and the isotopic composition of the incoming fluid.

4.5 Implication for mantle wedge hydration

Infiltration of the slab-derived fluid into the mantle wedge is important for subduction zone settings because mantle minerals are strongly depleted in volatiles. At equilibrium, a free aqueous fluid is not stable in the mantle wedge at T<650C until a fully hydrated mineral assemblage has formed (i.e. serpentine, chlorite, talc, and amphibole) (Manning, 2004). As shown above, the PT conditions of H2O release from the subducting slab, as well as the volume and the δ18O of the liberated H2O, can vary according to the geometry of the subduction zone and the composition of the subducting lithosphere (e.g. Hacker, 2008; Hernández-Uribe and Palin, 2019; Poli and Schmidt, 2002). The program PTLoop allows the PT conditions, amount and oxygen isotope composition of the released fluid to be calculated for the system of interest. The presented model uses an intermediate PT gradient that stabilizes lawsonite and results in a first significant fluid pulse at 65–70 km depth and a second pulse at ∼80 km depth. Below that depth, phengite is the main carrier of H2O. This implies that the majority of fluid is released in the forearc region, in agreement with previous studies investigating the dehydration of basaltic and sedimentary components of the slab (e.g. Baxter and Caddick, 2013; Kerrick and Connolly, 2001; Schmidt and Poli, 1994).

The influence of the slab-derived fluid on (1) the degree of hydration and (2) the δ18O modification of the overlying mantle rocks was estimated on the basis of the results presented above. A slab composed by the metabasalt(a) and the metasediment(t) (left column in Fig. 1a, assuming the real thickness of the slab to be 3 times the modelled one) subducting at 1 cm yr−1 was considered (Fig. 9a). This subduction rate represents a conservative estimate, considering that previous averages have proposed 4–5 cm yr−1 (Stern, 2002). Any faster subducting rate may shorten the timescale of the processes discussed below, allowing larger fluid fluxes in the mantle wedge over the same interval. Mechanical decoupling between the slab and the wedge and a steady-state cold mantle wedge are assumed (e.g. Abers et al., 2006; Hirauchi and Katayama, 2013; Wada et al., 2008). The fluid released by the slab at 500–520 C with a characteristic δ18O= 15.0 ‰ (Fig. 6a) was allowed to infiltrate into an initially dry peridotite (composition KLB-1 from Walter, 1998, simplified to the Fe–Mg–Al–Si system; Table S2.4 in Supplement S2) with T=570C and an initial δ18O= 5.5 ‰ (Eiler et al., 1997; Mattey et al., 1994). This simplified model ignores the dynamics of fluid infiltration and assumes pervasive flow. Given the column length of 1 m, a subduction rate of 1 cm yr−1 implies that, in 100 years, any fixed point (i.e. fixed PT conditions) at the slab–mantle interface receives the total amount of fluid that a single column can liberate in those conditions. Hence, in this example 4892.6 kg of water per 100 years (i.e. the amount released by the considered column at the chosen conditions, see above) infiltrate the mantle wedge. The released fluid first interacts with a small volume at the slab–mantle interface (V1 1000×1000×1 m; Fig. 9b) and once it has equilibrated and saturated V1, it infiltrates volume V2 (3000×2000×1 m; Fig. 9b). Both the volumes were scaled to 1:3 for the calculation in order to maintain the volume proportion with the modelled slab. The slab-derived fluid, which is in continuous supply during subduction of new material, infiltrates V1 and causes a progressive change in mineralogy from olivine + orthopyroxene +  garnet to serpentine + chlorite + minor olivine until the rock reaches saturation after 0.35 Myr of subduction. At this stage, V1 has a bulk δ18O of ∼8 ‰ that is significantly higher with respect to the initial mantle value (Fig. 9c). With ongoing subduction, the continuing release of water from new slab material under a static mantle drives the δ18O of the volume of mantle wedge toward higher values. The water that will infiltrate V2 has a δ18O that depends on (1) the δ18O of the slab-derived water and (2) the buffering capacity of V1. The same changes in mineral assemblage described for V1 occur also in V2, while the bulk δ18O of V2 increases more moderately than the one of V1 and reaches a δ18O of ∼6 ‰ after 0.75 Myr of subduction (Supplement S4).

Figure 9Case model for mantle wedge hydration. (a) Sketch of a subduction zone (crustal thickness and serpentinized mantle wedge dimensions are from Bostock et al., 2002). The subducting slab is composed of metabasalt(a) and metasediment(t) as shown in Fig. 1a, left column. (b) Geometry of the model. The blue arrow represents slab dehydration at 500–520 C. Abbreviations: V1 represents the volume of mantle rocks at the interface and V2 the surrounding volume (see text for details). (c) Plot of the bulk δ18O variations of V1 and V2 as a consequence of continuous slab dehydration over 1 Myr.


In the proposed model, most of the fluid is released by the slab at forearc depths. However, in most subduction zones no melting appears to occur in the forearc region and the serpentinite acts as the effective H2O absorber (Iwamori, 1998), recording the possible variation in δ18O induced by the slab-derived fluid. Progressive oceanward migration of the slab (“slab rollback”) has been regarded as an important mechanism acting in most active subduction zones (e.g. Heuret and Lallemand, 2005; Nakakuki and Mura, 2013). The rollback of the slab results in a lateral extension of the serpentinized wedge. As a consequence, the melt ascending below the arc can interact with serpentinized, high-δ18O mantle portions that were originally part of the forearc mantle and modify its original isotopic composition. High-δ18O arc lavas and melt in arc-setting peridotites have been described (e.g. phenocrysts in lavas from central Kamchatka: olivine δ18O= 5.8 ‰–7.1 ‰ and clinopyroxenes δ18O= 6.2 ‰–7.5 ‰, Dorendorf et al., 2000; New Guinea: silicate glass inclusions in olivine δ18O= 8.8 ‰–12.2 ‰, clinopyroxenes in metasomatized lehrzolite δ18O= 6.2 ‰–10.3 ‰, Eiler et al., 1998), but the mechanism of crustal contamination is still debated. Our results support the model proposed by Auer et al. (2009) that relates such high-δ18O lavas to the interaction between primitive basaltic melts with an uppermost mantle that was hydrated and enriched as part of the forearc mantle prior to trench migration.

4.6 Model applications and future directions

The presented approach has a broad range of applications for modelling fluid–rock interaction in different tectonic settings. We have presented here an example of subducted crust but the same principles apply also for regional metamorphism or hydrothermal systems. The model also provides new ways to quantify the degree of interaction of an external fluid within the same rock unit. We have shown that the observed effect on the δ18O of a rock of channelled vs. pervasive hydration is strictly coupled with the composition of the fluid source. Nevertheless, important insights can be given by linking observations of ideal cases with modelling even if the composition of the infiltrating fluid is not known a priori. For instance, the oxygen isotope composition of a fluid source can be retrieved when a variation in δ18O is observed within the same rock type from the more hydrated to the less hydrated portions, even in the absence of a clear presence of a vein or vein system. Alternatively, the degree of equilibration of the host rock around a vein can be calculated when the isotopic composition of the fluid is known. Therefore, the combination of mineral-scale in situ oxygen isotope analyses with major and trace element mapping will provide much more detailed information on quantitative element mobility during fluid–rock interaction.

Fluids play a central role as a catalyst for chemical reactions in rocks. Generally re-equilibration reactions occur only in the presence of fluids that either derive from breakdown of hydrous phases or from external sources (e.g. Airaghi et al., 2017; Cartwright and Barnicoat, 2003; Engi et al., 2018; Konrad-Schmolke et al., 2011; Rubatto and Angiboust, 2015). The program PTLoop calculates at which PT conditions the breakdown of hydrous phases occurs, and, consequently, metamorphic reactions and free fluid are expected. If fluid-driven reactions occurred in a rock at PT conditions where no release of internal fluid is expected, the role of an external fluid should be considered, and its amount and isotopic composition can be retrieved using the approach outlined in this paper.

Fluid–rock reactions in subducted lithosphere are likely to involve more complex open-system processes than the dehydration and rehydration considered for this model, driving silicate and carbonate dissolution, transport and re-precipitation (e.g. Ague and Nicolescu, 2013; Piccoli et al., 2016). Carbon release via the dissolution of calcium carbonate has been recognized to have important implications for CO2 release from subduction zones, and it is controlled by H2O-rich fluid infiltration (e.g. Ague and Nicolescu, 2013; Frezzotti et al., 2011; Gorman et al., 2006; Kerrick and Connolly, 2001). Future models may account for such variations in reactive major element bulk composition of the rocks along the PT evolution as a consequence of mineral net transfer reactions occurring simultaneously with fluid–rock exchange (Baumgartner and Valley, 2001), in addition to water liberation and mineral fractionation (see above). This would require additional considerations on the changes in the fluid isotopic composition during transfer of solute species through the fluid. The implementation of other fluid species beside H2O, such as CO2, could be assessed, provided that (1) reliable constraints on the oxygen isotope fractionation between these species and water or minerals are determined and (2)  their consistency with other available data is established. Other species such as CH4 and H2 do not contain any oxygen, and thus they are likely to be less relevant to this model.

5 Conclusions

We developed a user-friendly tool that combines equilibrium thermodynamic with oxygen isotope fractionation modelling for investigating the interaction between fluids and minerals in rocks during their metamorphic evolution. The program simulates along any given PT path the stable mineral assemblages, bulk δ18O and δ18O of stable phases and the amount and oxygen isotope composition of the fluid released.

The capabilities of the program PTLoop are illustrated by an application to subduction zones, but the presented modelling strategy can be applied to various metamorphic and tectonic settings. In this study, the chosen system represents a section of subducting oceanic crust composed by a lower layer of metabasalt and an upper layer of metasediments of carbonaceous or pelitic composition. The calculation follows a step-wise procedure along the chosen PT path. During the prograde evolution, any mineral and excess fluid can be fractionated from the reactive bulk composition.

Mineral fractionation and/or excess fluid loss produce only minor (i.e. < 1.0 ‰) shifts in the bulk δ18O of any lithology. Hence, the bulk δ18O of a rock that experienced a succession of such processes without interaction with external fluids is likely to be representative of its protolith composition. Variations in δ18O of stable phases due to mineral fractionation and/or excess fluid loss are also negligible (i.e. < 0.5 ‰), while the effect of temperature variation over a range of ∼150C on the mineral δ18O is phase dependent and may be significant (> 1.0 ‰).

Interaction with an external fluid of different oxygen isotope composition leads to shifts in bulk and mineral δ18O values according to the degree of fluid–rock interaction and δ18O difference between the rock and the fluid. Extremely large variations in bulk δ18O of ∼12 ‰ are calculated for the carbonate metasediment equilibrating with a fluid derived from a metabasalt with an initial hydrated MORB composition, while small variations of ∼3 ‰ are calculated for the terrigenous metasediment equilibrating with a fluid from a metabasalt that derives from an altered oceanic crust. When 50 % or more of the fluid deriving from dehydration of the metabasalts equilibrates with any of the overlying metasediments, the final δ18O of the fluid released by the system has a dominant sedimentary signature, with values between 12 ‰ and 18 ‰. Such fluids have δ18O values significantly higher than the mantle value (5.5 ‰) and have a great potential to modify the oxygen isotope composition of the mantle wedge at the slab–mantle interface. Extensive serpentinization and a δ18O increase of ∼2.5 ‰ are modelled at the interface already after 0.35 Myr of ongoing subduction.

PTLoop provides a powerful way to evaluate the effect of closed-system vs. open-system behaviour with respect to oxygen isotopes during the evolution of the rocks. Different degrees of interaction between the external fluids and the sink lithology can be simulated and the effects of internally vs. externally buffered fluids on the mineral paragenesis and on the mineral isotopic composition investigated.

Measured oxygen isotope compositions in minerals, intra-grain or bulk δ18O variations at different scales can be compared with the results of the model for specific scenarios. If the measured isotopic compositions are not consistent with the behaviour of a closed system, the presented approach can be used to determine feasible external fluid sources, to estimate the degree of fluid–rock interaction and the metamorphic conditions at which this happened. This modelling strategy can also assist in retrieving the oxygen isotope composition of a fluid source when a variation in δ18O is observed within the same rock type from the more hydrated to the less hydrated portions, even in the absence of a clear presence of a vein or vein system. Our model thus opens new avenues for mapping fluid pathways related to external fluid infiltration during the metamorphic evolution of the crust, with important consequences for element recycling in subduction zones and the investigation of fluid-induced earthquakes.

Code availability

A compiled version of the program PTLoop, the oxygen isotope fractionation database and the thermodynamic database used for this study are available at (Lanari and Vho, 2020).


The supplement related to this article is available online at:

Author contributions

AV developed the model and performed the calculations. PL supervised the software development and contributed to the code implementation. DR and JH contributed to formulate and design the model and to the interpretation of the results. AV prepared the paper with contributions from all co-authors. DR conceived the project and secured funding.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “Exploring new frontiers in fluids processes in subduction zones”. It is a result of the EGU Galileo conference “Exploring new frontiers in fluids processes in subduction zones”, Leibnitz, Austria, 24–29 June 2018.


We thank the conveners and the participants of the EGU Galileo conference “Exploring new frontiers in fluid processes in subduction zones” for constructive discussions. We also thank Ulrich Linden for the technical support and the reviewers Ralf Halama and Alberto Vitale Bovarone for the constructive comments.

Financial support

This research has been supported by the Swiss National Science Foundation to Daniela Rubatto (grant no. 200021_166280) and Jörg Hermann (grant no. 200021_169062).

Review statement

This paper was edited by Nadia Malaspina and reviewed by Ralf Halama and Alberto Vitale Brovarone.


Abers, G. A., van Keken, P. E., Kneller, E. A., Ferris, A., and Stachnik, J. C.: The thermal structure of subduction zones constrained by seismic imaging: Implications for slab dehydration and wedge flow, Earth Planet. Sc. Lett., 241, 387–397,, 2006. 

Ague J. J. and Nicolescu, S.: Carbon dioxide released from subduction zones by fluid-mediated reactions, Nat. Geosci., 7, 355–360,, 2014. 

Airaghi, L., Lanari, P., de Sigoyer, J., and Guillot, S.: Microstructural vs compositional preservation and pseudomorphic replacement of muscovite in deformed metapelites from the Longmen Shan (Sichuan, China), Lithos, 282, 262–280,, 2017. 

Albarède, F.: The survival of mantle geochemical heterogeneities, Washington DC American Geophysical Union Geophysical Monograph Series, 160, 27–46,, 2005. 

Alt, J. C., Muehlenbachs, K., and Honnorez, J.: An oxygen isotopic profile through the upper kilometer of the oceanic crust, DSDP Hole 504B, Earth Planet. Sc. Lett., 80, 217–229,, 1986. 

Auer, S., Bindeman, I., Wallace, P., Ponomareva, V., and Portnyagin, M.: The origin of hydrous, high-δ18O voluminous volcanism: diverse oxygen isotope values and high magmatic water contents within the volcanic record of Klyuchevskoy volcano, Kamchatka, Russia, Contrib. Mineral. Petr., 157, 209–230,, 2009. 

Barnes, J. D. and Straub, S. M.: Chorine stable isotope variations in Izu Bonin tephra: implications for serpentinite subduction, Chem. Geol., 272, 62–74,, 2010. 

Barnicoat, A. C. and Cartwright, I.: Focused fluid flow during subduction: oxygen isotope data from high-pressure ophiolites of the western Alps, Earth Planet. Sc. Lett., 132, 53–61,, 1995. 

Baumgartner, L. P. and Valley, J. W.: Stable isotope transport and contact metamorphic fluid flow, Rev. Mineral. Geochem., 43, 415–467,, 2001. 

Baxter, E. F. and Caddick, M. J.: Garnet growth as a proxy for progressive subduction zone dehydration, Geology, 41, 643–646,, 2013. 

Bebout, G. E. and Penniston-Dorland, S. C.: Fluid and mass transfer at subduction interfaces – The field metamorphic record, Lithos, 240, 228–258,, 2016. 

Bebout, G. E., Agard, P., Kobayashi, K., Moriguti, T., and Nakamura, E.: Devolatilization history and trace element mobility in deeply subducted sedimentary rocks: Evidence from Western Alps HP/UHP suites, Chem. Geol., 342, 1–20,, 2013. 

Bostock, M. G., Hyndman, R. D., Rondenay, S., and Peacock, S. M.: An inverted continental Moho and serpentinization of the forearc mantle, Nature, 417, 536–538,, 2002. 

Bowman, J. R., Willett, S. D., and Cook, S. J.: Oxygen isotopic transport and exchange during fluid flow: One-dimensional models and applications, Am. J. Sci., 294, 1–55,, 1994. 

Bulle, F., Rubatto, D., Ruggieri, G., Luisier, C., Villa, I. M., and Baumgartner, L.: Episodic hydrothermal alteration recorded by microscale oxygen isotope analysis of white mica in the Larderello-Travale Geothermal Field, Italy, Chem. Geol., 532, 119288,, 2019. 

Caddick, M. J., Konopásek, J., and Thompson, A. B.: Preservation of garnet growth zoning and the duration of prograde metamorphism, J. Petrol., 51, 2327–2347,, 2010. 

Cartwright, I. and Barnicoat, A. C.: Stable isotope geochemistry of Alpine ophiolites: a window to ocean-floor hydrothermal alteration and constraints on fluid–rock interaction during high-pressure metamorphism, Int. J. Earth Sci., 88, 219–235,, 1999. 

Cartwright, I. and Barnicoat, A. C.: Geochemical and stable isotope resetting in shear zones from Täschalp: constraints on fluid flow during exhumation in the Western Alps, J. Metamorph. Geol., 21, 143–161,, 2003. 

Cook-Kollars, J., Bebout, G. E., Collins, N. C., Angiboust, S., and Agard, P.: Subduction zone metamorphic pathway for deep carbon cycling: I. Evidence from HP/UHP metasedimentary rocks, Italian Alps, Chem. Geol., 386, 31–48,, 2014. 

de Capitani, C. and Brown, T. H.: The computation of chemical equilibrium in complex systems containing non-ideal solutions, Geochim. Cosmochim. Ac., 51, 2639–2652,, 1987. 

de Capitani, C. and Petrakakis, K.: The computation of equilibrium assemblage diagrams with Theriak/Domino software, Am. Mineral., 95, 1006–1016,, 2010. 

Diener, J. F. A., Powell, R., White, R. W., and Holland, T. J. B.: A new thermodynamic model for clino-and orthoamphiboles in the system Na2O–CaO–FeO–MgO–Al2O3SiO2H2O–O, J. Metamorph. Geol., 25, 631–656,, 2007. 

Dorendorf, F., Wiechert, U., and Wörner, G.: Hydrated sub-arc mantle: a source for the Kluchevskoy volcano, Kamchatka/Russia, Earth Planet. Sc. Lett., 175, 69–86,, 2000. 

Dubacq, B., Vidal, O., and De Andrade, V.: Dehydration of dioctahedral aluminous phyllosilicates: thermodynamic modelling and implications for thermobarometric estimates, Contrib. Mineral. Petr., 159, 159–174,, 2010. 

Eiler, J. M.: Oxygen isotope variations of basaltic lavas and upper mantle rocks, Rev. Mineral. Geochem., 43, 319–364,, 2001. 

Eiler, J. M., Farley, K. A., Valley, J. W., Hauri, E., Craig, H., Hart, S. R., and Stolper, E. M.: Oxygen isotope variations in ocean island basalt phenocrysts, Geochim. Cosmochim. Ac., 61, 2281–2293,, 1997. 

Eiler, J. M., McInnes, B., Valley, J. W., Graham, C. M., and Stolper, E. M.: Oxygen isotope evidence for slab-derived fluids in the sub-arc mantle, Nature, 393, 777–781,, 1998. 

Engi, M., Giuntoli, F., Lanari, P., Burn, M., Kunz, B., and Bouvier, A.-S.: Pervasive Eclogitization Due to Brittle Deformation and Rehydration of Subducted Basement: Effects on Continental Recycling?, Geochem. Geophy. Geosy., 19, 865–881,, 2018. 

Errico, J. C., Barnes, J. D., Strickland, A., and Valley, J. W.: Oxygen isotope zoning in garnets from Franciscan eclogite blocks: evidence for rock–buffered fluid interaction in the mantle wedge, Contrib. Mineral. Petr., 166, 1161–1176,, 2013. 

Evans, T. P.: A method for calculating effective bulk composition modification due to crystal fractionation in garnet-bearing schist: implications for isopleth thermobarometry, J. Metamorph. Geol., 22, 547–557,, 2004. 

Ferry, J. M. and Spear, F. S.: Experimental calibration of the partitioning of Fe and Mg between biotite and garnet, Contrib. Mineral. Petr., 66, 113–117,, 1978. 

Ferry, J. M., Kitajima, K., Strickland, A., and Valley, J. W.: Ion microprobe survey of the grain-scale oxygen isotope geochemistry of minerals in metamorphic rocks, Geochim. Cosmochim. Ac., 144, 403–433,, 2014. 

Frey, M., de Capitani, C., and Liou, J. G.: A new petrogenetic grid for low-grade metabasites, J. Metamorph. Geol., 9, 497–509,, 1991. 

Frezzotti, M. L., Selverstone, J., Sharp, Z. D., and Compagnoni, R.: Carbonate dissolution during subduction revealed by diamond-bearing rocks from the Alps, Nat. Geosci., 4, 703–706,, 2011. 

Friedman, I. and O'Neil, J. R.: Compilation of stable isotope fractionation factors of geochemical interest, in: Data of Geochemistry, 6th Edn., Geological Survey Professional Paper 440 KK,, 1977. 

Früh-Green, G. L., Scambelluri, M., and Vallis, F.: O–H isotope ratios of high pressure ultramafic rocks: implications for fluid sources and mobility in the subducted hydrous mantle, Contrib. Mineral. Petr., 141, 145–159,, 2001. 

Gale, A., Dalton, C. A., Langmuir, C. H., Su, Y., and Schilling, J.-G.: The mean composition of ocean ridge basalts, Geochem. Geophy. Geosy., 14, 489–518,, 2013. 

Gauthiez-Putallaz, L., Rubatto, D., and Hermann, J.: Dating prograde fluid pulses during subduction by in situ U–Pb and oxygen isotope analysis, Contrib. Mineral. Petr., 171, 1–20,, 2016. 

Gerdes, M. L., Baumgartner, L. P., Person, M., and Rumble, D.: One-and two-dimensional models of fluid flow and stable isotope exchange at an outcrop in the Adamello contact aureole, Southern Alps, Italy, Am. Mineral., 80, 1004–1019,, 1995a. 

Gerdes, M. L., Baumgartner, L. P., and Person, M.: Stochastic permeability models of fluid flow during contact metamorphism, Geology, 23, 945–948, 1995b. 

Gerya, T. V., Stöckhert, B., and Perchuk, A. L.: Exhumation of high-pressure metamorphic rocks in a subduction channel: A numerical simulation, Tectonics, 21, 6–1,, 2002. 

Giuntoli, F., Lanari, P., and Engi, M.: Deeply subducted continental fragments – Part 1: Fracturing, dissolution–precipitation, and diffusion processes recorded by garnet textures of the central Sesia Zone (western Italian Alps), Solid Earth, 9, 167–189,, 2018. 

Gorman, P. J., Kerrick, D. M., and Connolly, J. A. D.: Modeling open system metamorphic decarbonation of subducting slabs, Geochem. Geophy. Geosy., 7, 1–21,, 2006. 

Gregory, R. T. and Taylor Jr., H. P.: An oxygen isotope profile in a section of Cretaceous oceanic crust, Samail ophiolite, Oman: Evidence for δ18O buffering of the oceans by deep (> 5 km) seawater-hydrothermal circulation at mid-ocean ridges, J. Geophys. Res.-Sol. Ea., 86, 2737–2755,, 1981. 

Hacker, B. R.: H2O subduction beyond arcs, Geochem. Geophy. Geosy., 9, 1–24,, 2008. 

Hermann, J., Müntener, O., and Scambelluri, M.: The importance of serpentinite mylonites for subduction and exhumation of oceanic crust, Tectonophysics, 327, 225–238,, 2000. 

Hernández-Uribe, D. and Palin, R. M.: A revised petrological model for subducted oceanic crust: insights from phase equilibrium modelling, J. Metamorph. Geol., 37, 475–768,, 2019. 

Heuret, A. and Lallemand, S.: Plate motions, slab dynamics and back-arc deformation, Phys. Earth Planet. Int., 149, 31–51,, 2005. 

Higashino, F., Rubatto, D., Kawakami, T., Bouvier, A.-S., and Baumgartner, L. P.: Oxygen isotope speedometry in granulite facies garnet recording fluid/melt–rock interaction (Sør Rondane Mountains, East Antarctica), J. Metamorph. Geol., 37, 1037–1048,, 2018. 

Hirauchi, K. and Katayama, I.: Rheological contrast between serpentine species and implications for slab–mantle wedge decoupling, Tectonophysics, 608, 545–551,, 2013. 

Hoefs, J.: Stable isotope geochemistry, 7th Edn., Springer International Publishing, Switzerland, 389 pp.,, 2015. 

Holland, T. and Powell, R.: Thermodynamics of order-disorder in minerals; II, Symmetric formalism applied to solid solutions, Am. Mineral., 81, 1425–1437,, 1996. 

Holland, T. and Powell, R.: Activity–composition relations for phases in petrological calculations: an asymmetric multicomponent formulation, Contrib. Mineral. Petr., 145, 492–501,, 2003. 

Holland, T., Baker, J., and Powell, R.: Mixing properties and activity-composition relationships of chlorites in the system MgO-FeO-Al2O3-SiO2-H2O, Eur. J. Mineral., 10, 395–406,, 1998. 

Holland, T. J. B. and Powell, R.: An internally consistent thermodynamic data set for phases of petrological interest, J. Metamorph. Geol., 16, 309–343,, 1998. 

Iwamori, H.: Transportation of H2O and melting in subduction zones, Earth Planet. Sc. Lett., 160, 65–80,, 1998. 

John, T., Scambelluri, M., Frische, M., Barnes, J. D., and Bach, W.: Dehydration of subducting serpentinite: implications for halogen mobility in subduction zones and the deep halogen cycle, Earth Planet. Sc. Lett., 308, 65–76,, 2011. 

Kerrick, D. M. and Connolly, J. A. D.: Metamorphic devolatilization of subducted marine sediments and the transport of volatiles into the Earth's mantle, Nature, 411, 293–296,, 2001. 

Kohn, M. J.: Modeling of prograde mineral δ18O changes in metamorphic systems, Contrib. Mineral. Petr., 113, 249–261,, 1993. 

Konrad-Schmolke, M., O'Brien, P. J., de Capitani, C., and Carswell, D. A.: Garnet growth at high-and ultra-high pressure conditions and the effect of element fractionation on mineral modes and composition, Lithos, 103, 309–332,, 2008. 

Konrad-Schmolke, M., Zack, T., O'Brien, P., and Barth, M.: Fluid migration above a subducted slab – Thermodynamic and trace element modelling of fluid–rock interaction in partially overprinted eclogite-facies rocks (Sesia Zone, Western Alps), Earth Planet. Sc. Lett., 311, 287–298,, 2011. 

Kuhn, B. K., Reusser, E., Powell, R., and Günther, D.: Metamorphic evolution of calc-schists in the Central Alps, Switzerland, Schweiz. Miner. Petrog., 85, 175–190,, 2005. 

Lanari, P. and Duesterhoeft, E.: Modeling Metamorphic Rocks Using Equilibrium Thermodynamics and Internally Consistent Databases: Past Achievements, Problems and Perspectives, J. Petrol., 60, 19–56,, 2019. 

Lanari, P. and Engi, M.: Local bulk composition effects on metamorphic mineral assemblages, Rev. Mineral. Geochem., 83, 55–102,, 2017. 

Lanari, P. and Vho, A.: Fractionation of oxygen isotopes between minerals: a powerful tool for thermometry and investigating fluid-rock interaction, available at:, last access: January 2020. 

Lanari, P., Giuntoli, F., Loury, C., Burn, M., and Engi, M.: An inverse modeling approach to obtain P–T conditions of metamorphic stages involving garnet growth and resorption, Eur. J. Mineral., 29, 181–199,, 2017. 

Laurent, V., Lanari, P., Naïr, I., Augier, R., Lahfid, A., and Jolivet, L.: Exhumation of eclogite and blueschist (Cyclades, Greece): Pressure–temperature evolution determined by thermobarometry and garnet equilibrium modelling, J. Metamorph. Geol., 36, 769–798,, 2018. 

Liou, J. G., Tsujimori, T., Zhang, R. Y., Katayama, I., and Maruyama, S.: Global UHP metamorphism and continental subduction/collision: the Himalayan model, Int. Geol. Rev., 46, 1–27,, 2004. 

Manning, C. E.: The chemistry of subduction-zone fluids, Earth Planet. Sc. Lett., 223, 1–16,, 2004. 

Martin, L. A., Rubatto, D., Crépisson, C., Hermann, J., Putlitz, B., and Vitale-Brovarone, A.: Garnet oxygen analysis by SHRIMP-SI: Matrix corrections and application to high-pressure metasomatic rocks from Alpine Corsica, Chem. Geol., 374, 25–36,, 2014. 

Mattey, D., Lowry, D., and Macpherson, C.: Oxygen isotope composition of mantle peridotite, Earth Planet. Sc. Lett., 128, 231–241,, 1994. 

Mével, C.: Serpentinization of abyssal peridotites at mid-ocean ridges, C. R. Geosci., 335, 825–852,, 2003. 

Miller, J. A. and Cartwright, I.: Distinguishing between seafloor alteration and fluid flow during subduction using stable isotope geochemistry: examples from Tethyan ophiolites in the Western Alps, J. Metamorph. Geol., 18, 467–482,, 2000. 

Miller, J. A., Cartwright, I., Buick, I. S., and Barnicoat, A. C.: An O-isotope profile through the HP–LT Corsican ophiolite, France and its implications for fluid flow during subduction, Chem. Geol., 178, 43–69,, 2001. 

Muehlenbachs, K. and Clayton, R. N.: Oxygen isotope studies of fresh and weathered submarine basalts, Can. J. Earth Sci., 9, 172–184,, 1972. 

Nakakuki, T. and Mura, E.: Dynamics of slab rollback and induced back-arc basin formation, Earth Planet. Sc. Lett., 361, 287–297,, 2013. 

O'Neil, J. R. and Adami, L. H.: Oxygen isotope partition function ratio of water and the structure of liquid water, J. Phys. Chem., 73, 1553–1558,, 1969. 

O'Neil, J. R. and Taylor, H. P.: The oxygen isotope and cation exchange chemistry of feldspars, Am. Mineral., 52, 1414–1437, 1967. 

Padrón-Navarta, J. A., Hermann, J., Garrido, C. J., Sánchez-Vizcaíno, V. L., and Gómez-Pugnaire, M. T.: An experimental investigation of antigorite dehydration in natural silica-enriched serpentinite, Contrib. Mineral. Petr., 159, 25–42,, 2010. 

Padrón-Navarta, J. A., Sánchez-Vizcaíno, V. L., Hermann, J., Connolly, J. A., Garrido, C. J., Gómez-Pugnaire, M. T., and Marchesi, C.: Tschermak's substitution in antigorite and consequences for phase relations and water liberation in high-grade serpentinites, Lithos, 178, 186–196,, 2013. 

Page, F. Z., Kita, N. T., and Valley, J. W.: Ion microprobe analysis of oxygen isotopes in garnets of complex chemistry, Chem. Geol., 270, 9–19,, 2010. 

Page, F. Z., Essene, E. J., Mukasa, S. B., and Valley, J. W.: A garnet–zircon oxygen isotope record of subduction and exhumation fluids from the Franciscan Complex, California, J. Petrol., 55, 103–131,, 2013. 

Page, F. Z., Cameron, E. M., Flood, C. M., Dobbins, J. W., Spicuzza, M. J., Kitajima, K., Strickland, A., Ushikubo, T., Mattinson, C. G., and Valley, J. W.: Extreme oxygen isotope zoning in garnet and zircon from a metachert block in mélange reveals metasomatism at the peak of subduction metamorphism, Geology, 47, 655–658,, 2019. 

Parra, T., Vidal, O., and Agard, P.: A thermodynamic model for Fe–Mg dioctahedral K white micas using data from phase-equilibrium experiments and natural pelitic assemblages, Contrib. Mineral. Petr., 143, 706–732,, 2002. 

Peacock, S. A.: Fluid processes in subduction zones, Science, 248, 329–337,, 1990. 

Peacock, S. M.: The importance of blueschist eclogite dehydration reactions in subducting oceanic crust, Geol. Soc. Am. Bull., 105, 684–694, 1993. 

Penniston-Dorland, S. C., Kohn, M. J., and Manning, C. E.: The global range of subduction zone thermal structures from exhumed blueschists and eclogites: Rocks are hotter than models, Earth Planet. Sc. Lett., 428, 243–254,, 2015. 

Piccoli, F., Vitale Brovarone, A., Beyssac, O., Martinez, I., Ague, J. J., and Chaduteau, C.: Carbonation by fluid–rock interactions at high-pressure conditions: implications for carbon cycling in subduction zones, Earth Planet. Sc. Lett., 445, 146–159,, 2016. 

Plank, T.: The chemical composition of subducting sediments, Treatise on Geochemistry, 4, 607–629,, 2014. 

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

Poli, S. and Schmidt, M. W.: Petrology of subducted slabs, Annu. Rev. Earth Pl. Sc., 30, 207–235,, 2002. 

Powell, R. and Holland, T. J. B.: On thermobarometry, J. Metamorph. Geol., 26, 155–179,, 2008. 

Putlitz, B., Matthews, A. M., and Valley, J. W.: Oxygen and hydrogen isotope study of high-pressure metagabbros and metabasalts (Cyclades, Greece): implications for the subduction of oceanic crust, Contrib. Mineral. Petr., 138, 114–126,, 2000. 

Rubatto, D. and Angiboust, S.: Oxygen isotope record of oceanic and high-pressure metasomatism: a P–T–time–fluid path for the Monviso eclogites (Italy), Contrib. Mineral. Petr., 170, 1–16,, 2015. 

Rumble, D.: Stable isotope fractionation during metamorphic devolatilization reactions, Rev. Mineral. Geochem., 10, 327–353, 1982. 

Russell, A. K., Kitajima, K., Strickland, A., Medaris, L. G., Schulze, D. J., and Valley, J. W.: Eclogite-facies fluid infiltration: constraints from δ18O zoning in garnet, Contrib. Mineral. Petr., 165, 103–116,, 2013. 

Scambelluri, M. and Philippot, P.: Deep fluids in subduction zones, Lithos, 55, 213–227,, 2001. 

Scambelluri, M., Fiebig, J., Malaspina, N., Müntener, O., and Pettke, T.: Serpentinite subduction: implications for fluid processes and trace-element recycling, Int. Geol. Rev., 46, 595–613,, 2004. 

Schmidt, M. W. and Poli, S.: The stability of lawsonite and zoisite at high pressures: Experiments in CASH to 92 kbar and implications for the presence of hydrous phases in subducted lithosphere, Earth Planet. Sc. Lett., 124, 105–118,, 1994. 

Sharp, Z.: Principles of stable isotope geochemistry, 2nd Edn., University of New Mexico, 384 pp.,, 2017. 

Sharp, Z. D. and Barnes, J. D.: Water-soluble chlorides in massive seafloor serpentinites: a source of chloride in subduction zones, Earth Planet. Sc. Lett., 226, 243–254,, 2004. 

Siron, G., Baumgartner, L., Bouvier, A.-S., Putlitz, B., and Vennemann, T.: Biotite reference materials for secondary ion mass spectrometry 18O∕16O measurements, Geostand. Geoanal. Res., 41, 243–253,, 2017. 

Spandler, C. and Hermann, J.: High-pressure veins in eclogite from New Caledonia; Implications for fluid migration and seismic activity in subduction zones, Geochim. Cosmochim. Ac., Supplement, 69, A664,, 2005. 

Spear, F. S.: Metamorphic fractional crystallization and internal metasomatism by diffusional homogenization of zoned garnets, Contrib. Mineral. Petr., 99, 507–517,, 1988. 

Spear, F. S., Pattison, D. R. M., and Cheney, J. T.: The metamorphosis of metamorphic petrology, in The Web of Geological Sciences: Advances, Impacts, and Interactions II, Geol. Soc. Am., 523, 31–73,, 2017. 

Staudigel, H.: Chemical fluxes from hydrothermal alteration of the oceanic crust, in: Treatise on Geochemistry, 2nd Edn., edited by: Turekian, H. D. H. K., Elsevier, Oxford, UK, 583–606, 2014. 

Staudigel, H., Plank, T., White, B., and Schmincke, H.-U.: Geochemical fluxes during seafloor alteration of the basaltic upper oceanic crust: DSDP Sites 417 and 418, Subduction: top to bottom, 96, 19–38,, 1996.  

Stern, R. J.: Subduction zones, Rev. Geophys., 40, 1–40,, 2002. 

Sun, S.-S. and McDonough, W. F.: Chemical and isotopic systematics of oceanic basalts: implications for mantle composition and processes, Geol. Soc., London, Special Publications, 42, 313–345,, 1989. 

Syracuse, E. M., van Keken, P. E., and Abers, G. A.: The global range of subduction zone thermal models, Phys. Earth Planet. Int., 183, 73–90,, 2010. 

Tracy, R. J.: Compositional zoning and inclusions in metamorphic minerals. Characterization of Metamorphism through Mineral Equilibria, Rev. Mineral., 10, 355–397, 1982. 

Valley, J.: Stable isotope geochemistry of metamorphic rocks, Rev. Mineral. Geochem., 16, 445–489, 1986. 

Vho, A., Lanari, P., and Rubatto, D.: An internally-consistent database for oxygen isotope fractionation between minerals, J. Petrol.,, 2020. 

Vidal, O., Lanari, P., Munoz, M., Bourdelle, F., de Andrade, V., and Walker, J.: Deciphering temperature, pressure and oxygen-activity conditions of chlorite formation, Clay Miner., 51, 615–633,, 2016. 

Vielzeuf, D., Veschambre, M., and Brunet, F.: Oxygen isotope heterogeneities and diffusion profile in composite metamorphic-magmatic garnets from the Pyrenees, Am. Mineral., 90, 463–472,, 2005a. 

Vielzeuf, D., Champenois, M., Valley, J. W., Brunet, F., and Devidal, J. L.: SIMS analyses of oxygen isotopes: matrix effects in Fe–Mg–Ca garnets, Chem. Geol., 223, 208–226,, 2005b. 

Wada, I., Wang, K., He, J., and Hyndman, R. D.: Weakening of the subduction interface and its effects on surface heat flow, slab dehydration, and mantle wedge serpentinization, J. Geophys. Res.-Sol. Ea., 113, B04402,, 2008. 

Walter, M. J.: Melting of garnet peridotite and the origin of komatiite and depleted lithosphere, J. Petrol., 39, 29–60,, 1998. 

White, W. M. and Klein, E. M.: Composition of the Oceanic Crust, Treatise on Geochemistry, 2nd Edn., 4, 457–496,, 2014. 

Whitney, D. L. and Evans, B. W.: Abbreviations for names of rock-forming minerals, Am. Mineral., 95, 185–187,, 2010. 

Williams, M. J.: Tracing Fluids from Seafloor to Deep Subduction Environments: An In-situ Geochemical Investigation of Fluid-Mobile Elements in Oceanic Crust, Ph.D thesis, Australian National University, Canberra,, 325 pp., 2019. 

Zack, T. and John, T.: An evaluation of reactive fluid flow and trace element mobility in subducting slabs, Chem. Geol., 239, 199–216,, 2007. 

Zheng, Y. F.: Oxygen isotope fractionation in metal monoxides, Mineral. Mag., 58, 1000–1001,, 1994. 

Short summary
This study presents an approach that combines equilibrium thermodynamic modelling with oxygen isotope fractionation modelling for investigating fluid–rock interaction in metamorphic systems. An application to subduction zones shows that chemical and isotopic zoning in minerals can be used to determine feasible fluid sources and the conditions of interaction. Slab-derived fluids can cause oxygen isotope variations in the mantle wedge that may result in anomalous isotopic signatures of arc lavas.