This contribution presents an approach and a computer program (GrtMod) for numerical simulation of garnet evolution based on compositions of successive growth zones in natural samples. For each garnet growth stage, a new local effective bulk composition is optimized, allowing for resorption and/or fractionation of previously crystallized garnet. The successive minimizations are performed using the Nelder–Mead algorithm; a heuristic search method. An automated strategy including two optimization stages and one refinement stage is described and tested. This program is used to calculate pressure–temperature (P–T) conditions of crystal growth as archived in garnet from the Sesia Zone (Western Alps). The compositional variability of successive growth zones is characterized using standardized X-ray maps and the program XMapTools. The model suggests that Permian garnet cores crystallized under granulite-facies conditions at T > 800 °C and P = 6 kbar. During Alpine times, a first garnet rim grew at eclogite-facies conditions (650 °C, 16 kbar) at the expense of the garnet core. A second rim was added at lower P (∼11 kbar) and 630 °C. In total, garnet resorption is modeled to amount to ∼9 vol% during the Alpine evolution; this value is supported by our observations in X-ray compositional maps.

1. Introduction

An inverse approach to a scientific problem involves the determination of the causal factors that satisfy a set of observations. In metamorphic rocks of a given bulk-rock composition (CBR), pressure (P) and temperature (T) conditions determine the stable mineral assemblage, thus they constitute the causal factors. In this case, observations are the coexisting phases defining the mineral assemblage, their compositions and volume proportions. However, Gibbs free energy minimization, the method classically used to model such an assemblage, is a forward technique (e.g., de Capitani & Petrakakis, 2010). Resulting equilibrium phase diagrams are strictly based on assemblages predicted by Gibbs free energy minimization for a given CBR, i.e., the composition of a rock volume devoid of compositional heterogeneities. Such diagrams combined with mineral isopleths have been intensively used to estimate P–T conditions in metamorphic rocks by comparing results of the model with observations.

Porphyroblasts are large crystals surrounded by a matrix of finer-grained minerals and are of central interest because they often preserve a chemical and textural record of metamorphic processes and conditions. For instance, garnet porphyroblast in low- to moderate-grade metamorphic rocks often are compositionally zoned, and in favorable cases such zoning can be used to infer the part of the P–T history during which garnet grew in a sample (e.g., Spear & Selverstone 1983; Spear et al., 1984, 1991a; Konrad-Schmolke et al., 2005; Gaidies et al., 2006, 2008b; Cheng & Cao, 2015). This method works for relatively large garnet crystals provided that metamorphic duration is not unusually long and that thermal maximum reached by the sample was below ∼700 °C. In such cases, intracrystalline diffusion is slow enough (Yardley, 1977; Caddick et al., 2010; Stowell et al., 2011) to preserve fine compositional differences, and the compositional zoning in garnet is likely to reflect growth conditions. However, variations in P–T conditions and isolation of early garnet growth zones imply a gradual change in CBR of the reactive part of the rock (Evans, 2004). In cases where a significant amount of compositionally zoned garnet is preserved (>2–3 vol%), a single equilibrium phase diagram cannot be used to retrieve successive P–T stages. During the last 25 years four main approaches were developed to overcome this problem, all of which are essentially based on thermodynamic equilibrium theory:

  • – DiffGibbs program (Spear et al., 1991b) allows for prediction of the chemical-zoning pattern of garnet. It accounts for intragranular diffusion in garnet operating simultaneously with net-transfer and exchange reactions during garnet growth. Along a specific P–T trajectory, garnet is assumed to be in equilibrium with a given set of phases. The limit of this approach is that it does not test the stability of the assemblage for the given CBR.

  • – The second approach consists of successive forward models for which CBR is manually altered to account for material sequestered in garnet cores, as analyzed and mapped by electron microprobe (Marmo et al., 2002; Tinkham & Ghent, 2005; Caddick et al., 2007).

  • – Forward modeling of garnet zoning and coexisting phases (including mineral modes) is realized for an arbitrary P–T path using Gibbs free energy minimization, eventually comparing at each stage the predicted composition of garnet with the observed zoning (e.g., Konrad-Schmolke et al., 2008; Robyr et al., 2014). Garnet compositions and volumes produced at each step are fractionated from the CBR, thus providing a new effective bulk composition. That composition refers to the volume domain in the rock over which thermodynamic equilibrium is established during one increment of garnet growth. The first limitation of such models is the arbitrary choice of the P–T trajectory. To improve the results of such models, Moynihan & Pattison (2013) provided an automated inverse strategy to derive the “best” P–T trajectory by minimizing a misfit parameter, basically the weighted differences between measured and model compositional profiles. Once a P–T point is found, a search begins for the next P–T point, with a model composition that best matches the next point on the garnet transect. Similarly, Vrijmoed & Hacker (2014) proposed a brute-force computational method (inverse technique as well) to determine the best P–T trajectory by minimizing the differences between predicted and measured garnet compositional profiles along different trajectories. The fundamental limit of both Moynihan & Pattison (2013) and Vrijmoed & Hacker (2014) approaches is that garnet growth is assumed to occur continuously, and no garnet resorption is taken into account. In reality, part of the fractionating garnet may continue to react during the next P–T point and to being dissolved. Evidence of garnet resorption is commonly visible in compositional maps (see, e.g., Figs. 3a,c and 4a in Moynihan & Pattison, 2013). Again, this may produce significant changes in effective composition, depending on the amount resorbed, which is typically not visible in a sample.

  • – The program Theria_G (Gaidies et al., 2008a) allows for numerical simulation of porphyroblast nucleation and garnet growth in a given volume of rock for any defined P–T–time (P–T–t) path. Two major modules are used by Theria_G, (1) the Gibbs free energy minimization routine of Theriak (de Capitani & Brown, 1987) and (2) a model describing intragranular multi-component diffusion. In contrast to the previous techniques, Theria_G simulates the formation of an entire population of garnet with variable grain size using a forward model. The arbitrary choice of the P–T trajectory is again a severe limitation of this model. However, Moynihan & Pattison (2013) used the approach described above to derive the best P–T trajectory, which is subsequently defined as input in Theria_G models.

Addressing some of the limitations of existing approaches, this study proposes an alternative strategy and a computer program, GrtMod (available at http://grtmod.petrochronology.org), to model garnet growth during successive P–T stages based on natural compositional records. To improve the control data, garnet compositions of successive growth zones are characterized from standardized X-ray maps (see details in Lanari et al., 2014). The P–T conditions as well as proportions of garnet resorption are optimized by GrtMod at each step to match the modeled and measured compositions. The model presented in this report differs from those reviewed above in that it is strongly based on the observation of preserved garnet growth zones in natural rocks; no provision for intracrystalline diffusion is made. GrtMod is written in MATLAB© and interacts with Theriak (de Capitani & Brown, 1987) using the extension Theriak_D (Duesterhoeft & de Capitani, 2013).

2. GRTMOD strategy

The strategy behind GrtMod consists of using unconstrained nonlinear optimization to find the minimum of an objective function in n-dimensional space (N ≥ 2). The variables to be optimized are P, T, and, for stages Si (i > 1), the volume fractions of all previous garnet growth zones that are fractionated from the CBR. These volume fractions may decrease in the course of modeling because of garnet resorption, which possibly affects the earlier growth zones to variable extents. The objective function used reflects the deviation of the model composition from the measured compositions. This approach critically relies on the characterization of the composition of successive garnet growth zones.

2.1. Growth stages of garnet and corresponding variables

As discussed in the introduction, in low- to moderate-grade metamorphic rocks (T < 700 °C), for relatively large garnet crystals and assuming that metamorphic duration is not unusually long, intracrystalline diffusion is slow and the successive garnet compositions are likely to reflect changes in equilibrium conditions only. Basically, the absence of zoning caused by intracrystalline diffusion can be verified by sharp compositional boundaries between the successive growth zones (Figs. 1 and 2).

Fig. 1

Compositional maps of garnet in a polymetamorphic eclogitic micaschist from the Sesia Zone. Data obtained using XMapTools: maps of end-member proportions generated using the external function Gar-StructForm: (a) grossular; (b) almandine, (c) pyrope; (d) spessartine. (e) Maps of the compositional groups generated using the module Chem2d (Lanari et al., 2014). Domains used to derive the average composition of Grt1, Grt2 and Grt3 are outlined by dashed line. (Online version in color.)

Fig. 1

Compositional maps of garnet in a polymetamorphic eclogitic micaschist from the Sesia Zone. Data obtained using XMapTools: maps of end-member proportions generated using the external function Gar-StructForm: (a) grossular; (b) almandine, (c) pyrope; (d) spessartine. (e) Maps of the compositional groups generated using the module Chem2d (Lanari et al., 2014). Domains used to derive the average composition of Grt1, Grt2 and Grt3 are outlined by dashed line. (Online version in color.)

Fig. 2

Zoning profile of garnet end-member proportions along AB transect (see Fig. 1): (a) grossular; (b) almandine; (c) pyrope; note different scales. Fractures in garnet core (Grt1) are indicated by white arrows in (a); the corresponding compositions are neglected in further calculations (black circles). Mean compositions of each garnet group are reported (solid black lines) together with standard deviation (1σ, dashed black line). Black bands show mineral inclusions in garnet.

Fig. 2

Zoning profile of garnet end-member proportions along AB transect (see Fig. 1): (a) grossular; (b) almandine; (c) pyrope; note different scales. Fractures in garnet core (Grt1) are indicated by white arrows in (a); the corresponding compositions are neglected in further calculations (black circles). Mean compositions of each garnet group are reported (solid black lines) together with standard deviation (1σ, dashed black line). Black bands show mineral inclusions in garnet.

In the present study, the growth history of garnet is divided into discrete events, defined as “growth stages” (Si, i = 1,…,n, see Table 1). An individual growth stage is a short event occurring at given Pi and Ti during which a garnet volume fraction vi grows with a homogeneous composition Cgrti.

For any growth stage Si, the variables to be optimized are its specific P and T (Pi and Ti) as well as the volume fraction of previous garnet to be fractionated from the CBR (vi,j= 1,…, vi,j=i−1). The number of variables at stage i thus is:  
Ni=2+(i1).
(1)
For the first stage, the CBR of the sample is assumed to be equal to the composition of the reactive part of the system; hence CBR serves as input for the forward models. For subsequent stages (i > 1), the local effective bulk composition, CLEBi, is calculated as follows:  
CLEBi=CBRj=1i=1vi,jρgrtiρrocki1Cgrti(1j=1i=1vi,j)ρmtxi1ρrocki1,
(2)
where Cgrti is the measured composition (in oxide wt%) of the garnet growth zone i; ρgrti is its model density, and ρrocki1 and ρmtxi1 are the average rock density and matrix density (i.e., all phases except garnet) from the preceding stage; vi,j is the volume fraction of garnet crystallized during stage j that is fractionated from the CBR at stage i. As some garnet of stage i is preserved in the present-day sample despite possible resorption, the following condition is applied by GrtMod at each growth zone:  
0<vi,jvi1,j,
(3)
The resorption of garnet j at stage i (i > j) is expressed in vol% of garnet in the bulk rock and is given by:  
ri,j=(vi1,jvi,j)×100.
(4)
The total resorption of garnet rgrt is defined as  
rgrt=ijri,j.
(5)
In order to compute the local effective bulk composition from garnet volume fractions – to fractionate the previous garnet from the CBR – the rock density of the previous stage is required (Eq. (2)). The density of the rock at stage i−1 (ρrocki1)is calculated assuming zero porosity:  
ρrocki1=j=1i1vi,jρgrti+ρmtxi1(1j=1i1vi,j).
(6)

The volume fractions used in this study are expressed as fractions of the entire system, i.e., the rock sample. Consequently the volume of garnet predicted by GrtMod to be stable for stages i > 1 must be corrected for the size of the subsystem being considered at each stage by the model (in Fig. 3, this corresponds to the blue, green, and red domains for stages 1, 2, and 3, respectively).

Fig. 3

Sketches illustrating garnet growth during three hypothetical stages that include resorption prior to growth stages 2 and 3. The colored domains are the subsystems being considered at each stage by the model. For stages 2 and 3, two sketches are displayed; the first one (left) shows the amount of resorption of previously formed garnet that is used to estimate the local effective bulk of that specific stage. The second sketch (right) shows the growth of a new garnet rim modeled using Gibbs free energy minimization. (Online version in color.)

Fig. 3

Sketches illustrating garnet growth during three hypothetical stages that include resorption prior to growth stages 2 and 3. The colored domains are the subsystems being considered at each stage by the model. For stages 2 and 3, two sketches are displayed; the first one (left) shows the amount of resorption of previously formed garnet that is used to estimate the local effective bulk of that specific stage. The second sketch (right) shows the growth of a new garnet rim modeled using Gibbs free energy minimization. (Online version in color.)

2.2. Objective function

The objective function is composed of a loss function generating the number L0 and a cost function generating the number C0. Both parameters are used to quantify the amount by which the predicted garnet composition deviates from the measured values. The loss function is defined as:  
L0(Pi,Ti,vi,j=1,,vi,j=i1)=k=1m(fkmeasuredfkmodel)2ωk,
(7)
where fkmeasured and fkmodel are the proportions of end-member k, as measured or modeled, and ωk is the weighting factor of the corresponding element. The end-member proportions of grossular (fGrs), pyrope (fPrp), almandine (fAlm), and spessartine (fSps) are calculated from Ca, Mg, Fe, and Mn abundances, expressed in number of atoms per formula unit (apfu). However, the proportions fGrs, fPrp, fAlm and fSps are not known with the same precision, hence a weighting factor is required. It takes into account the relative analytical uncertainty of each end-member proportion by its variance σk2 (σk: standard deviation), hence the weighting factor is defined as  
ωk=1σk2.
(8)
In mapping conditions, the precision (1σ) of the electron microprobe measurement of one pixel composition is estimated using a Poisson law (e.g., Lanari et al., 2014) 
p=1n,
(9)
where n is the number of photons reaching the detector during a single measurement. The number of recorded counts n is corrected for dead-time bias of the detector. For multiple independent measurements fk of the same X-ray flux f, the variance is close to  
σk2=n.
(10)
This obvious relationship is derived from equation (9) and may be tested by calculating the variance of single measurements of the same composition. For a large set of measurements of a homogeneous material – such as a set of pixels of X-ray maps – the precision calculated from the variance matches the precision of the single-pixel estimate made using equation (9). By combining equations (8) and (10), it is possible to estimate the weighting factor of an end-member k using the number of recorded counts nk of the corresponding element using the relationship  
ωk=1nk.
(11)
Consequently, the loss function used in this study is  
L0(Pi,Ti,vi,j=1,,vi,j=i1)=k=1m(fkmeasuredfkmodel)2nk.
(12)
The number L0 generated by the loss function is minimized to derive the best set of variables (maximum likelihood solution). However, this number is not intuitively representing the deviation between model and measured compositions because of the weighting factor. This is the reason why a cost function is also part of the objective function. In contrast to the loss function, the cost function generates the number C0, which is intuitively representing the quality of the solution. The cost function is defined as:  
C0(Pi,Ti,vi,j=1,,vi,j=i1)=k=1m(fkmeasuredfkmodel)2,
(13)

C0 is the least square of the deviations between the model and measured end-member proportions without taking into account the uncertainty on the measurement. For garnet with four end-members, a value of C0 < 0.04 indicates a good fit of the model. For garnet showing low Mn-content (<1 wt %), only three end-members may be used to describe its composition. In such a case, a threshold value of 0.03 is set for C0 (see application example below).

2.3. Minimization procedure

The problem addressed in this study is a nonlinear optimization problem for which the derivatives are not known. Consequently a heuristic search method has to be used; the Nelder–Mead technique (Nelder & Mead, 1965), implemented in the MATLAB© function fminsearch, was selected. It is critical for the user to understand the technique of minimization in order to evaluate the limits of this approach. A complete method description is available in Nelder & Mead (1965) or when using the help function in MATLAB©. This method uses the concept of simplex, a polytope of n + 1 vertices in n dimensions. The n + 1 values of the objective function L at the vertices are ordered and the position of the centroid is calculated of all n points, except for the worst point n + 1. Then the algorithm computes the values of the objective function at the reflected point, the expanded point and the contracted point of the worst point. If one of the previous values is smaller, the corresponding point replaces the worst point and the optimization continues, else a reduction step of the simplex is done. This procedure is repeated until convergence to a local minimum. The best solution found is a local minimum and may differ significantly from the global minimum of the objective function, which represents the maximum likelihood solution. Consequently the optimization must be done using different initial guesses. An attempt to provide an automated procedure is proposed in the next section. Advantages and disadvantages of this automated procedure are pointed out in the discussion.

2.4. Toward an automated optimization strategy

The optimization strategy depends on the number of variables being optimized. For instance, for the first stage S1, only two variables (P1, T1) are optimized, and the procedure is straightforward. For all subsequent stages (Si>1), the problem is more complex because of the additional compositional variables (Fig. 3). The procedure is divided into three phases: optimization1, optimization2, and auto-refinement. During optimization1, successive P–T minimizations are carried out from different starting guesses in order to determine the global minimum within the P–T window (Fig. 4a). Optimization2 refines the garnet fraction variables (vi,j= 1→i−1), as well as P–T (Fig. 4b). The starting guess for optimization2 is the best P–T couple obtained during optimization1. A go fast mode is available to begin directly optimization2 from user’s favorite P–T initial guess. Finally, the auto-refinement phase checks the local variability of the cost function (C0 value) in order to provide a relative uncertainty on the P–T estimate (Fig. 5). A complete description of these three phases is presented in the Appendix 1.

Fig. 4

Sketches illustrating the automated procedure. (a) Pressure– temperature diagram within TminTmax, PminPmax. The G1,ri are the four initial guesses (r = 1:4). (b) Pressure (P)–temperature (T)–composition (X) diagram describing the two-stage optimization. Optimization1 is carried out in P–T space (red), optimization2 in P–T–v space (blue). In this example C0 of S1,2i is smaller than C0 of S1,1i and is selected as the best P–T couple obtained during optimization1. (Online version in color.)

Fig. 4

Sketches illustrating the automated procedure. (a) Pressure– temperature diagram within TminTmax, PminPmax. The G1,ri are the four initial guesses (r = 1:4). (b) Pressure (P)–temperature (T)–composition (X) diagram describing the two-stage optimization. Optimization1 is carried out in P–T space (red), optimization2 in P–T–v space (blue). In this example C0 of S1,2i is smaller than C0 of S1,1i and is selected as the best P–T couple obtained during optimization1. (Online version in color.)

Fig. 5

Sketch illustrating the refinement procedure, which checks the local variability of the cost function. Two hypothetic solutions 1 and 2 are found for stage i. D1(n1), D2(n2), …, D8(n8) give the directions (D) and the number of steps calculated along them (n). dP and dT are the P and T increments. The gray areas correspond to domains with C0 < TC0. For a given direction D, the refinement stops if C0 > TC0. (Online version in color.)

Fig. 5

Sketch illustrating the refinement procedure, which checks the local variability of the cost function. Two hypothetic solutions 1 and 2 are found for stage i. D1(n1), D2(n2), …, D8(n8) give the directions (D) and the number of steps calculated along them (n). dP and dT are the P and T increments. The gray areas correspond to domains with C0 < TC0. For a given direction D, the refinement stops if C0 > TC0. (Online version in color.)

2.5. Compositional characterization of growth zones

As discussed above, during a single growth stage the fractionation of the effective composition caused by garnet growth as well as small changes in P and T conditions are neglected. These two assumptions are critical and demand extensive characterization of the compositional variations of the studied sample. It is strongly recommended to use high-resolution quantitative compositional maps of garnet end-members to define the successive growth zones (see, e.g., areas used in Fig. 1). The quantitative X-ray mapping technique is useful to measure the compositional variability of metamorphic minerals at the thin-section scale (e.g., Lanari et al., 2012, 2013). Local compositional variations caused by surface crystal defects, late re-equilibration with inclusions or fractures are not taken into account; any analyses showing mixed analyses close to the contact between two growth zones are discarded. In Fig. 2, for example, the intermediate compositions between the defined growth stages are not being considered. It is not excluded that they may result from protracted growth with a minor change in P–T conditions or kinetic-controlled growth. This approach critically relies on the chemical information stored in the natural sample and, therefore, on the quality of the microprobe measurements. No a priori assumption is made on the P–T conditions of each stage. However, some garnet compositions from some growth stages may have been totally resorbed during later stages. If the composition is not preserved in the present-day specimen, it is obvious that P–T conditions cannot be retrieved using this approach.

3. Benchmarking test for a sample showing typical prograde garnet zoning

To benchmark GrtMod, P–T conditions of a garnet in the San Emigdio Schist (sample 06SE23 from Chapman et al., 2011) were estimated in a system simplified to SiO2-Al2O3-FeO-MnO-MgO-CaO-Na2O-K2O-H2O using the same thermodynamic dataset (TC321p2.txt) and CBR as in that study. This example was selected because (i) the authors demonstrated that the effects of intracrystalline diffusion were limited to narrow zones (∼10 μm; Fig. 7 in Chapman et al., 2011), and (ii) there is no clear evidence of resorption (see Fig. 6a,b in Chapman et al., 2011). Thus the automated optimization strategy of GrtMod is expected to predict incremental growth of garnet along the prograde P–T history without resorption.

The zoning profile reported in Fig. 5b of Chapman et al. (2011) was divided into four successive growth zones: Grt1 (0.37 vol%, Alm36Prp2Grs26Sps36); Grt2 (0.12 vol%, Alm57Prp3Grs33Sps7); Grt3 (4.2 vol%, Alm67Prp7Grs25Sps1); and Grt4 (2.31 vol%, Alm64Prp9Grs26Sps1). The volume fraction of each growth zone was calculated assuming a total of 7 vol% of garnet being produced during the prograde P–T history (value taken from Fig. 11b of Chapman et al., 2011). As explained above (see Sect. 2), the procedure requires a subdivision into discrete stages of growth, and four stages where selected based on the zoning profile.

The benchmark results are reported in Fig. 6. The model predicts that garnet grows along a prograde trajectory with increasing T and P from ∼450 °C to ∼630 °C and 4 kbar to 11 kbar. It is important to point out here that the model does not predict any resorption of the previous growth zones (Fig. 6b), in line with the conclusions of Chapman et al. (2011), and the model results were obtained without any intervention from the user. The predicted zoning profiles are compared with the measured one in Fig. 6c,d. Although Grt2 is slightly underestimated in the model (0.03 vol% instead of 0.12 vol%), the predicted profile shape perfectly matches the observations. The residuals are very low (C0 between 0.007 and 0.025) resulting in an excellent match of the modeled compositions.

Fig. 6

Results of the benchmark test (see text for details). (a) P–T diagram with the solutions obtained for the four stages; best solutions are highlighted. Reported in gray is the P–T path obtained by Chapman et al. (2011). (b) Volume fraction of garnet predicted by the best model for stages 1, 2, 3, and 4. (c) Zoning profiles (thick lines) and predictions from this study (thin) for almandine and pyrope (left); grossular and spessartine (right). (Online version in color.)

Fig. 6

Results of the benchmark test (see text for details). (a) P–T diagram with the solutions obtained for the four stages; best solutions are highlighted. Reported in gray is the P–T path obtained by Chapman et al. (2011). (b) Volume fraction of garnet predicted by the best model for stages 1, 2, 3, and 4. (c) Zoning profiles (thick lines) and predictions from this study (thin) for almandine and pyrope (left); grossular and spessartine (right). (Online version in color.)

4. Sample description and compositional mapping

The studied sample FG12-157 is an eclogitic garnet-bearing micaschist from the Sesia Zone in the Italian Western Alps (see Supplementary Material S6 for photographs). It was selected from a collection of ∼10 samples showing similar garnet resorption features because it illustrates well the strengths and weaknesses of the automated approach. Other samples, some of which show more P–T stages or no resorption (Giuntoli, 2016), will be presented in a subsequent study. Sample FG12-157 was collected at Lillianes in the Lys Valley in Italy (X = 409683; Y = 5054033 ED 1950 UTM Zone 32N). This micaschist is made of quartz (40 vol%), phengite (30 vol%), garnet (15 vol%,) glaucophane (6 vol%), and epidote (4 vol%), with minor chlorite, albite, rutile, zircon, titanite, ilmenite, and graphite (all together of about 5 vol%). A strong eclogitic foliation is marked by phengite, glaucophane and allanite; it was subsequently deformed into open folds. Garnet-grain size ranges from 200 μm up to several mm. Microscopically, large garnet porphyroblasts systematically show a bright core surrounded by a dark rim. The small grains, however, are dark crystals with features similar to the rims of porphyroblasts. The cloudy appearance of the dark garnet is mostly due to fine rutile inclusions. Glaucophane and phengite inclusions are frequent in the dark rim, whereas only quartz is found in the bright core. In the matrix, glaucophane shows characteristic core to rim zoning, with more strongly pleochroic rims (darker blue compared to core) reflecting higher iron contents. Some glaucophane grains are rimmed by green amphibole, indicating local and limited retrogression, with minor chlorite and albite reflecting greenschist-facies conditions. Allanite shows REE-rich cores and locally a clinozoisite rim (10 μm). Other accessories are graphite, zircon, and rutile overgrown by titanite and followed by an ilmenite rim.

Electron probe microanalyses (EPMA) were performed using a JEOL JXA-8200 superprobe at the Institute of Geological Sciences (University of Bern). Following the procedure described in Lanari et al. (2013, 2014), the EPMA session is divided into two steps, i.e., the measurement of point analyses and X-ray compositional maps, both in wavelength-dispersive mode (WDS). Analytical conditions for point analyses were 15 kV accelerating voltage, 20 nA specimen current, and 40 s dwell time (including 2×10 s of background measurement). Nine oxide components were measured, using synthetic and natural standards: wollastonite (SiO2), anorthite (Al2O3, CaO), almandine (FeO), spinel (MgO), orthoclase (K2O), albite (Na2O), ilmenite (TiO2), and tephroite (MnO). Analytical conditions for X-ray maps were 15 kV accelerating voltage, 100 nA specimen current, and a dwell time of 200 ms/pixel. Nine elements (Si, Ti, Al, Fe, Mn, Mg, Na, Ca and K) were measured at the specific wavelength in two passes. Intensity maps were standardized using spot analyses as internal standard. X-ray maps were processed using XMapTools 2.2.1 (Lanari et al., 2014). The average composition of each growth zone was calculated from the map pixels selected (Fig. 1). Analytical uncertainties (derived using Eq. (9)) were propagated through the structural formulae computation using a Monte-Carlo simulation. Compositional data and corresponding analytical uncertainties are reported in Table 2. The sum of the analytical uncertainties on fPrp, fGrs and fAlm is about 0.03. This value is used to define the value of STOL.

5. Results

5.1. Garnet composition and texture

Garnet exhibits complex zoning as shown by compositional maps of end-member proportions (Fig. 1). The cores are Prp- and Alm-rich and Grs-poor (Alm68–70Prp24–26Grs3–5Sps1–3, Table 2). Cores have lobate edges suggesting resorption. A new garnet enriched in Grs and depleted in Prp and Alm (Alm66–70Prp17–21Grs10–14Sps1–3) fills up numerous fractures (Fig. 1). Two distinct over-growths surround the apparently porphyroclastic cores: a first rim (Alm63–65Prp18–20Grs15–17Sps1) and a second rim (Alm57–61Prp14–16Grs24–28Sps1) (Fig. 1e). Sample textures indicate that the second rim grew on internally and externally resorbed portions of the first rim, i.e., the second rim is observed directly surrounding the core as well as the first rim. This observation suggests that partial resorption of garnet core plus growth of the first rim may have occurred before or during growth of the second rim (third stage in Fig. 3). The first rim is identical in composition to garnet that seals hairline fractures in the core. This observation supports the sequence of crystallization proposed here. Resorption of garnet core is common in polymetamorphic rocks, sometimes leading to the formation of mushroom-shaped and atoll garnet (Robyr et al., 2014 and references therein). Based on these textural and compositional relationships, three growth stages are defined: stage 1 corresponding to the growth of garnet core (Grt1), stage 2 for the first rim (Grt2), and stage 3 for the second rim (Grt3).

5.2. Thermodynamic models

For this application example, the thermodynamic dataset of Berman (1988) with subsequent updates collected in JUN92.bs (distributed with Theriak-Domino 03.01.2012) was used together with the following solution models: Berman (1990) for garnet; Fuhrman and Lindsley (1988) for feldspar; Meyre et al. (1997) for omphacite; Keller et al. (2005) for white mica, and ideal mixing models for amphibole (Mäder & Berman, 1992; Mäder et al., 1994), epidote and chlorite (Hunziker, 2003). As the studied garnet contains <1 wt% MnO, restricted to parts of the core (suggesting heterogeneous distribution of the Mn-rich precursors), the Mn component is ignored and the system simplified to SiO2-TiO2-Al2O3-FeO-MgO-CaO-Na2O-K2O-H2O. For other cases, where Mn-rich garnet is modeled, the MnO component must be added to the system (e.g., benchmark test in Sect. 3). The anhydrous CBR, determined by XRF, comprises SiO2 (60.36 wt%), TiO2 (1.03 wt%), Al2O3 (16.51 wt%), FeO (7.95 wt%), MgO (3.29 wt%), CaO (2.09 wt%), Na2O (1.19 wt%), and K2O (3.57 wt%). Because of the lack of experimental data and suitable ferric end-members in solid-solution models, Fe3+ was ignored. All Gibbs free energy minimizations were carried out assuming a saturating pure H2O fluid. The amount of H2O predicted at high P is in line with the measured LOI (2.01 wt%) in the present-day sample. End-member mineral abbreviations used throughout text and figures are from Whitney & Evans (2010).

The model is restricted to a P–T window between 500–900 °C and 5–20 kbar. Minimum garnet abundances, i.e., those preserved in the present-day sample, were fixed at 4 vol% for Grt1, 3 vol% for Grt2, and 4 vol%for Grt3. Details regarding the results printed out by GrtMod and the modeled assemblages are reported in Supplementary Material S1–S5 (linked to this article and freely available online at the GSW website of the journal: http://eurjmin.geoscienceworld.org), corresponding to the stages described in the following subsections. Input and output values of selected variables are reported in Table 4.

5.2.1. Stage 1

During stage 1, only P1 and T1 are optimized. Four starting guesses G1,11, G1,21, G1,31, G1,41 were defined at 600 °C–8 kbar, 600 °C–16 kbar, 800 °C–8 kbar, and 800 °C–16 kbar, respectively, following the procedure described in Fig. 4. It is instructive to follow the iterations step by step and to describe the results provided in Table 3.

The first minimization ( G1,11 from 600 °C–8 kbar) converges to a minimum at 851 °C and 6.03 kbar for a C0 value of 0.021. A solution S1,11 is retained for the following because for this first minimization C0 < 0.03. Model fAlm and fGrs differ from the measured values by ∼0.01 each. 10.5 vol% of garnet is predicted to be stable at that stage. The second minimization ( G1,21 from 600 °C– 16 kbar) converges to a different minimum at 674 °C and 17.49 kbar with a C0 value of 0.107. As C0 is much higher, model fAlm and fGrs are quite different from the measured values (model, 0.62 and 0.11; measured, 0.70 and 0.04) and no solution is saved because C0 > 0.03. The third minimization ( G1,31 from 800 °C–8 kbar) converges to a minimum at 851 °C and 6.02 kbar for a C0 value of 0.021. This result is fairly similar to solution S1,11, indeed minimizations 1 and 3 converge to the same local minimum. The fourth minimization ( G1,41 from 800 °C– 16 kbar) converges to a minimum at 899 °C and 6.14 kbar for a C0 value of 0.015. A solution S1,21 is obtained (C0 < 0.03). 7.14 vol% of garnet is predicted to be stable at that stage. The second solution (S1,21) has a smaller value of C0(Co1,21<Co1,11) and is selected as the best solution for stage 1: Sbest1=S1,21 (see Fig. 7; Tables 3 and 4).

Fig. 7

Results for the application example (eclogitic micaschist from the Sesia Zone). (a) P–T diagram where all the solutions for the three stages are reported. The best solutions (see text) are highlighted. Note that the alternative solutions for stages 1 and 3 (not used in this model) are reported in gray. (b) Volume fraction of garnet predicted by the best model during stages 1, 2 and 3. (Online version in color.)

Fig. 7

Results for the application example (eclogitic micaschist from the Sesia Zone). (a) P–T diagram where all the solutions for the three stages are reported. The best solutions (see text) are highlighted. Note that the alternative solutions for stages 1 and 3 (not used in this model) are reported in gray. (b) Volume fraction of garnet predicted by the best model during stages 1, 2 and 3. (Online version in color.)

5.2.2. Stage 2 – automated procedure

For stage 2, the optimization is divided into two steps: optimization1 and optimization2. During optimization1, a fixed amount of garnet Grt1 (7.14 vol%), corresponding to the amount predicted during stage 1, is initially fractionated from the CBR. The P–T conditions of the four starting guesses G1,12, G1,22, G1,32, G1,42 are the same as for stage 1 (see above). The first minimization ( G1,12 from 600 °C–8 kbar) converges to a minimum at 719 °C and 9.66 kbar with C0 = 1.74×10−4. A temporary solution S1,12 is stored; 5.95 vol% of garnet is predicted to be stable at that stage. The second minimization ( G1,22 from 600 °C– 16 kbar) converges to a minimum at 647 °C and 16.03 kbar with C0 = 4.15×10−3. A temporary solution S1,22 is stored; 9.02 vol% of garnet is predicted to be stable at that stage. The third minimization ( G1,32 from 800 °C–8 kbar) converges to a solution S1,32 at 719 °C and 9.66 kbar with C0 = 1.25×10−4. This solution is close to S1,22 with a slightly smaller residual. The fourth minimization ( G1,42 from 800 °C–16 kbar) converges to a solution S1,422 similar to S1,22. Optimization1 of stage 2 shows that for the same CBR (i.e., without resorption) Grt2 is predicted stable at 719 °C and 9.66 kbar and 647 °C and 16.03 kbar. The automated algorithm selects S1,32 as the best solution based on the C0 values (Sbest2=S1,32).

The P–T conditions of initial guesses of optimization2 are fixed at 719 °C and 9.66 kbar. A new variable v2,1 corresponding to the quantity of garnet Grt1 crystallized during stage 1 and fractionated from the bulk composition during stage 2 is introduced in the variable list of the objective function. Three starting guesses are defined assuming no resorption (v2,1 = 7.14 vol%), strong resorption (v2,1 = 4 vol%) and moderate resorption (v2,1 = 5.57 vol%). The garnet volume fraction used as input for the second starting guess corresponds to the amount of garnet that is preserved in the present-day sample. The first minimization ( G2,12 of 7.14 vol%) converges to a minimum at 719 °C and 9.65 kbar and a final v2,1 of 7.13 vol%, with a C0 value of 6.55×10−5 and 5.95 vol% of Grt2. The second minimization ( G2,22 of 4 vol%) converges to a minimum at 694 °C and 9.49 kbar and a final v2,1 of 4.35 vol% with a C0 value of 4.44×10−5 and 9.44 vol% of Grt2. The third minimization ( G2,32 of 5.57 vol%) converges to a minimum at 707 °C and 9.57 kbar and a final v2,1 of 5.86 vol% with a C0 value of 2.46×10−5 and 7.53 vol% of Grt2. These results suggest that Grt2 is modeled between 694 °C and 719 °C and between 9.49 and 9.65 kbar. The increase in T is associated with a decrease in resorption (from 2.78 vol % to zero) and to a decrease in the fraction of Grt2 (from 9.44 to 5.95 vol%). In this case different scenarios can be selected for the second stage. However, as discussed below (see Sect. 6.3), this solution is not the most likely based on the mineral inclusions captured during the growth of this first rim.

5.2.3. Stage 2 – go fast mode

The results for stage 2 discussed above show that the automated procedure selects the best solution at the end of optimization1 and consequently can ignore a solution with different P–T values. The go fast mode is used with initial P–T being taken from S1,22 at 650 °C and 16 kbar. The user defines moderate resorption to stabilize garnet at such conditions. The minimization converges to a minimum at 647 °C and 15.95 kbar and a final v2,1 of 5.92 vol% with C0 = 8.12×10−6 and 10.40 vol% of Grt2. From a statistical point of view, this solution at higher P is better than those found with the automated procedure (see Sect. 5.2.2). However, both P–T conditions allow for precise modeling of the observed garnet compositions. In contrast to the low-P solutions, the P–T conditions obtained at high P with and without resorption are similar (ΔT = 0.37 °C, ΔP = 0.09 kbar). In this case, resorption of 1.21 vol% of Grt1 improves the quality of the model and is selected as the most likely solution for stage 2.

5.2.4. Stage 3 – automated procedure

For stage 3, the optimization is again divided into two steps: optimization1 and optimization2. During optimization1, a fixed quantity of garnet Grt1 and Grt2 (5.92 vol% and 6.00 vol%) is fractionated from the CBR. For Grt2 optimization1 is carried out assuming 4.40 vol% resorption, the minimum required to stabilize garnet with a composition similar to Grt3. The first minimization ( G1,13 of 600 °C–8 kbar) does not converge to a solution because the amount of garnet produced is too small (1.8 vol% predicted whereas the minimum amount of Grt3 is fixed at 4 vol% observed in the sample). The second minimization ( G1,23 from 600 °C–16 kbar) converges to a minimum at 637 °C and 13.33 kbar for C0 = 1.71×10−4 solution. A temporary S1,13 is stored; 6.17 vol% of garnet is predicted stable at that stage. The third minimization ( G1,33 from 800 °C–8 kbar) does not converge because it predicts only 2.2 vol% of garnet, far below the minimum amount for Grt3. The fourth minimization converges to the same minimum as S1,13 with a slightly smaller C0 value of 1.62×10−4. A temporary solution S1,23 is stored. The algorithm selects S1,23 as the best solution based on the C0 values (S1,best3=S1,23).

The P–T conditions of initial guesses of optimization2 are fixed at 637 °C and 13.33 kbar. Two new variables v3,1 and v3,2 are to be optimized; they correspond to the quantities of garnet Grt1 and Grt2 fractionated from the bulk composition during stage 3. Three starting guesses are defined assuming either no more resorption than given by the input value (v3,1 = 5.92 vol% and v3,2 = 6.00 vol%), strong resorption (v3,1 = 4.00 vol% and v3,2 = 3.00 vol%), and moderate resorption (v3,1 = 4.96 vol% and v3,2 = 4.50 vol%). In this case initial v3,2 is set at 6.00 vol% whereas the volume fraction of Grt2 at stage 2 is 10.40 vol%. As only a small volume fraction of Grt2 is preserved in the present-day sample (i.e., much less than was produced in stage 2), strong resorption is expected to occur during stage 3. The first minimization converges to a minimum (S2,13) at 637 °C and 13.33 kbar and v3,1 and v3,2 of 5.92 and 6.00 vol% with C0 = 9.72×10−5 and 6.16 vol% Grt3. The second minimization converges to a minimum (S2,23) at 632 °C and 11.01 kbar and v3,1 and v3,2 of 4.86 and 3.40 vol % with C0 = 2.68×10−6 and 8.42 vol% of Grt3. The third minimization converges to a minimum (S2,33) at 637 °C and 13.33 kbar and v3,1 and v3,2 of 5.16 and 4.62 vol% with C0 = 5.89×10−6 and 8.27 vol% Grt3. Two solutions are found at slightly different pressures (13.37 kbar and 11.01 kbar). The second solution S2,23 is selected here because it is considered more likely based on the C0 value and because it matches observed mineral proportions and textural observations better, suggesting stronger resorption of the first rim and core before crystallization of the second rim.

6. Discussion

6.1. Intergranular diffusion and global equilibrium within domains

In this study garnet growth is modeled based on CBR and assuming thermodynamic equilibrium to be achieved at the millimeter scale. Component transport in the rock matrix through an intergranular medium is assumed to be fast relative to garnet growth, minimizing chemical potential gradients in the matrix. For the studied sample, compositional maps show that garnet in quartz-rich layers recorded a zoning that is similar to that in phengite-rich layers. Such observations and the excellent match of model compositions support the assumption of global equilibrium through an intergranular medium during each individual growth stage. Characteristic diffusion distance for Al in an intergranular medium saturated with hydrous fluid at 650 °C is >1 cm for a time >1 Myr (Carlson, 2010), allowing for homogenization of the composition at the sample scale. However, it is well known that in some cases, porphyroblast growth can lead to development of local chemical heterogeneities generating changes in nutrient production rates (Carlson et al., 1995). In such cases, our model will not be suitable (Carlson et al., 2015), and a diffusion-controlled model should be used (see Konrad-Schmolke et al., 2005; Schwartz et al., 2011; Ketcham & Carlson 2012).

6.2. Heuristic search method and domains with local minimum

As this study deals with non-linear problems (see Sect. 2.3) requiring a heuristic search method, the Nelder–Mead technique was selected (Nelder & Mead, 1965). However, at the end of a single minimization it is not possible to ensure that the minimum found is the global minimum. One way to ensure that a convergence point is a global minimum is to map the objective-function with high P–T resolution. An algorithm to compute such P–T maps of C0 and L0 values for a given effective bulk composition is provided in GrtMod.

As an example, the C0 map in the P–T range 500–900 °C and 5–20 kbar was computed for stage 1 using the CBR (Fig. 8a,b). This map exhibits the shape of the cost function C for stage1 (seeSect.5.2.1). Two distinct regions with local minima are found by the automated function: the first region at high P (HP; 674 °C, 17.49 kbar, C0 = 0.107) and the second region at high T (HT; 850–900 °C, ∼6 kbar, C0 < 0.03). Every single minimization starting on one side of the ridge separating the two low regions (dashed line in Fig. 8a) converges to the HP domain. Those starting on the opposite flank of the ridge converge to the HT domain. For stage 1 the automated procedure finds the global minimum within the HT domain. This first example shows that the shape of such objective function can be complex. Hence it is crucial to run many successive minimizations from different starting guesses.

Fig. 8

2D and 3D P–T maps of the cost function C for stage 1 (a, b) and stage 2 without any resorption (c, d). The four starting guesses are reported in the 2D diagrams only. The dashed line outlines the ridge between the two low regions that show local minima (see text). (Online version in color.)

Fig. 8

2D and 3D P–T maps of the cost function C for stage 1 (a, b) and stage 2 without any resorption (c, d). The four starting guesses are reported in the 2D diagrams only. The dashed line outlines the ridge between the two low regions that show local minima (see text). (Online version in color.)

6.3. Automated strategy [1]: limitation of multiple minima and solution finding

Similarly, the C0 map in the range 500–900 °C and 5–20 kbar was computed for stage 2, assuming no resorption of Grt1 (Fig. 8c,d). The effective CBR is computed by subtracting from CBR the amount of Grt1 produced during stage 1. Two regions with local minima are found by the automated function during optimization1: one at HT (3 solutions around ∼720 °C and ∼6.6 kbar for a best C0 of 1.25×10−4), and the other at HP (647 °C and 16 kbar with C0 = 4.15×10−4). In that case, only the second minimization converges toward HP domain because the starting guess is located on the other flank of the ridge separating the two low regions. As the lower C0 without resorption is found within the HT domain, 647 °C and 16 kbar are selected as starting guess for optimization2. In such a case, the HP domain is not investigated during optimization2 (i.e., with resorption) by the automated function. However, this can be done using the go fast mode (see Sect. 5.2.3). The results described above demonstrate that for stage 2 the best solution from optimization2 (C0 = 4.15×10−6) is found in the HP domain. The automated function converges to a local minimum, which has distinct P–T conditions compared to those of the global minimum (HP). However, there is little difference in C0 between both solutions, and Grt2 composition can be accurately modeled at 720 °C/6.6 kbar and 647 °C/16 kbar.

Stage 2 shows that garnet alone can provide ambiguous results in the framework of deriving P–T conditions of one single metamorphic stage. In such cases, it is crucial to incorporate the study of the coexisting phases. For example, numerous inclusions of phengite in Grt2 suggest that it coexisted with Si-rich phengite (Si4+ = 3.33 ± 0.02 apfu, XMg = 0.76 ± 0.02) (Table 5). The K-rich white-mica composition predicted by the model at 647 °C/16 kbar is Si4+ = 3.34 and XMg = 0.79, whereas at 720 °C/6.6 kbar it is Si4+ = 3.17 and XMg = 0.68. As expected, the Si content in phengite increases with increasing P. Modeled phengite composition nicely matches the measured composition at 647 °C and 16 kbar. The study of coexisting phases, such as phengite in this example, strongly supports Grt2 growth during a HP stage.

6.4. Automated strategy [2]: P–T–X minimization

The two-step optimization proposed in this study is expected to work when the shape of the objective function does not change significantly with garnet resorption. The investigated sample exhibits evidence of strong garnet resorption (Fig. 1) and is thus well suited to demonstrate the effectiveness of this automated method. The overall goal of optimization1 is to find out the best P–T starting guess in a simple two-variable problem.

TheC0 maps in the P–Trange550–700 °C and 12–20 kbar were computed for stage 2 assuming 0%, 42% and 84% resorption of garnet core (Fig. 9). The P–T position of the best solution from optimization1 (black line in Fig. 9) is very similar to values that assume intermediate and strong resorption (red lines in Fig. 9). This example demonstrates that, for a restricted P–T range with a single minimum, the two-step optimization is an elegant strategy to solve the problem.

Fig. 9

3D P–T diagrams of the cost function C for stage 2 with variable resorption (a – 0%, b – 42%, c – 84% of garnet core). The black line marks the P–T position of the best solution from optimization1. The red lines mark the P–T position of the best solutions assuming moderate (42%) and strong (84%) resorption of garnet core. (Online version in color.)

Fig. 9

3D P–T diagrams of the cost function C for stage 2 with variable resorption (a – 0%, b – 42%, c – 84% of garnet core). The black line marks the P–T position of the best solution from optimization1. The red lines mark the P–T position of the best solutions assuming moderate (42%) and strong (84%) resorption of garnet core. (Online version in color.)

6.5. P–T stages recorded in garnet from a polymetamorphic micaschist

The model predicts that garnet core (Grt1) crystallized under granulite-facies conditions at T > 800 °C and ∼6 kbar. This result is in line with other estimates available for the same area (Lardeaux & Spalla, 1991; Rebay & Spalla, 2001; Giuntoli, 2016). Such HT/low-P metamorphic conditions are common in the area and were recorded during the Permian (Rebay & Spalla, 2001). The first rim (Grt2) is Alpine and grew under eclogite-facies conditions (650 ± 50 °C, 16 ± 2.5 kbar). Similar HP conditions have been proposed for nearby areas of the Sesia Zone (e.g., Konrad-Schmolke et al., 2011; Regis et al., 2014). The model predicts Grt2 to grow at the expense of Grt1. The shape of Grt1 remnants with lobate edges supports this result. However, it is not possible to establish precisely when resorption occurred. It happened after stage 1 and either before or during stage 2. The composition of Grt1 is very different from any observed in comparable rocks containing typical Alpine prograde garnet. The P–T conditions obtained for Grt1 suggest that it formed at granulite-facies conditions, prior to Alpine orogeny, most likely during the Permian. At these HT conditions, the protolith must have been largely dehydrated, hence a stage of rehydration must be invoked to explain the development of the mica-rich Alpine eclogite assemblage. A scenario of HP hydration that triggered the dissolution of Grt1 and the precipitation of Grt2 seems plausible, and it may explain why substantial reaction overstepping (prior to hydration) occurred (50–150 °C and 2–4 kbar, depending on the prograde trajectory). The second Alpine rim (Grt3) grew at lower P, estimated at 11 ± 2 kbar and 632 ± 50 °C. The two solutions reported in Fig. 7 show similar P–T ranges. This late stage may be associated with phengite rims and crossitic amphibole.

6.6. Importance of incorporating resorption in thermodynamic models

Our results from the textural analysis predict resorption of up to 36 vol% of the total garnet produced occurring at >11 kbar. During stage 3, the model implies that 70 vol% of Grt2 was resorbed. Such partial resorption has a strong effect on the effective CBR used in modeling. It is fair to ask what a classical model would predict, such as those discussed in the introduction (Gaidies et al., 2008a; Konrad-Schmolke et al., 2008; Moynihan & Pattison 2013; Vrijmoed & Hacker 2014), all of which do account for fractionation of garnet, but only for the amounts produced, not those resorbed. The performance of such models has been tested using three different model variants:

  • – Model T-1 (Fig. 10b): This model is computed used Theriak and is based on fractionation (100% of garnet produced) along a retrograde P–T path involving five steps between 15.95 kbar (647 °C) and 11.01 kbar (632 °C). These values were chosen based on thermobarometric results of this study (Fig. 10a) and on petrological evidence (e.g., phengite inclusions; see Sect. 6.3). The relict garnet core Grt1 was initially fractionated from the CBR in order to generate a suitable effective bulk composition for stage 2 (Fig. 10b). A limitation of this test is that the P–T trajectory was arbitrarily chosen.

    Fig. 10

    Zoning profiles (thick lines) and model compositions (thin) for almandine, pyrope and grossular contents of the Alpine rims. The models were generated using (a) GrtMod and the strategy described in this study, (b) Theriak along a fixed P–T loop from 15.95 kbar (647 °C) and 11.01 kbar (632 °C), and (c, d) GrtMod with a strategy similar to that of Moynihan & Pattison (2013) (MP). In MP-1, successive P–T optimization is done using the garnet composition of the next point on the zoning profile (black squares in c) without any consideration of the garnet volume produced. In MP-2, the successive P–T optimizations are done using the garnet composition of the point on the zoning profile that corresponds to the volume of garnet previously produced. For both cases ∼11 vol% of garnet is produced in three stages (labeled Grt2, Grt3 and Grt4). For each model, the P–T conditions are summarized in a small box. Zoning profiles of the Alpine rims in Fig. 2 are scaled and translated into fractional volume in order to represent the minimum volume fraction of Grt2 and Grt3 produced. This estimate was based on the compositional maps reported in Fig. 1 and on maps of other porphyroblasts from the same sample (Giuntoli, 2016).

    Fig. 10

    Zoning profiles (thick lines) and model compositions (thin) for almandine, pyrope and grossular contents of the Alpine rims. The models were generated using (a) GrtMod and the strategy described in this study, (b) Theriak along a fixed P–T loop from 15.95 kbar (647 °C) and 11.01 kbar (632 °C), and (c, d) GrtMod with a strategy similar to that of Moynihan & Pattison (2013) (MP). In MP-1, successive P–T optimization is done using the garnet composition of the next point on the zoning profile (black squares in c) without any consideration of the garnet volume produced. In MP-2, the successive P–T optimizations are done using the garnet composition of the point on the zoning profile that corresponds to the volume of garnet previously produced. For both cases ∼11 vol% of garnet is produced in three stages (labeled Grt2, Grt3 and Grt4). For each model, the P–T conditions are summarized in a small box. Zoning profiles of the Alpine rims in Fig. 2 are scaled and translated into fractional volume in order to represent the minimum volume fraction of Grt2 and Grt3 produced. This estimate was based on the compositional maps reported in Fig. 1 and on maps of other porphyroblasts from the same sample (Giuntoli, 2016).

  • – Model MP-1 (Fig. 10c): To avoid this arbitrariness, the strategy described by Moynihan & Pattison (2013) was used. For the first composition of the zoning profile the best conditions (P1T1) are found, then a second point is analyzed and finds P2T2, etc. For each step, the garnet model composition best matches the next point of the zoning profile. In MP-1, successive P–T optimizations are performed using the garnet composition of the subsequent point of the zoning profile (black squares in Fig. 10c); however, the amount (volume) of garnet produced is not considered.

  • – Model MP-2 (Fig. 10d): Same optimization as MP-1, but in this case the successive P–T optimizations are done using the garnet composition of the point on the zoning profile that corresponds to the previously produced volume of garnet (black squares in Fig. 10d).

Models MP-1 and MP-2 were computed using GrtMod with an option that prevents garnet resorption. In both cases ∼11 vol% of garnet is produced in three stages (labeled Grt2, Grt3 and Grt4 in Fig. 10c,d); this amount corresponds to the total volume of Alpine garnet found in the sample. Both models MP-1 and MP-2 retrieve the best P–T path for the given modeling conditions. In model T-1, garnet is no longer predicted stable after the second P–T step (for the local effective bulk composition).

The importance of considering resorption in forward thermodynamic models is evident when the results of the three classical models (with fractionation only) are compared with our reference model that includes resorption and fractionation. The reference model with two growth stages and partial resorption (GrtMod in Fig. 10a) matches the observed zoning profile as well as the volume fractions of each growth zone. All of the classical models (T-1,MP-1and MP-2) fail to reproduce the observed zoning profile. Furthermore, the fraction of Grt2 is always overestimated at the first stage of growth (∼16 kbar). For Grt3, the models MP-1 and MP-2 predict different P–T scenarios. MP-1 implies HP garnet because the model always tries to fit the first Alpine rim composition that is not satisfactory for the subsequent rims (Grt3 and Grt4 in Fig. 10c). In contrast, model MP-2 does match the composition of Grt3 but the corresponding volume fraction is seriously underestimated, whereas Grt2 is overestimated. These discrepancies are caused by the absence of resorption. Nevertheless, the P–T conditions predicted by MP-2 are similar to those of the reference model.

This simple example demonstrates how important it is to consider the volume of garnet, not only the shape of the compositional profile. By comparing predicted and observed volume fractions in the sample, the amount of growth and resorption can be estimated. For samples that experienced no resorption, the models tested here produce the same result (Fig. 6). Of course the comparison between samples and models is limited since the amount of garnet remaining from each growth zone after subsequent resorption is only a minimum of the garnet produced at that stage. Therefore, the success of a model covering all stages of growth and resorption needs to be judged by comparing all of their compositions and modal amounts.

7. Conclusions

In this study, we provide a strategy and a computer program, GrtMod, to model garnet growth through successive stages by minimizing differences between measured and modeled compositions as predicted for a given effective bulk composition. Gibbs free energy minimization is used to obtain the model results. During garnet growth, the previous growth zones can either be fractionated from the bulk rock composition or be partially resorbed.

The shape of the objective function may be complex, sometimes showing two distinct local minima at different P–T conditions (see Fig. 8). An automated strategy is proposed, but the results strongly rely on the first P–T optimization. The example with two solutions shows that it is crucial to compare the measured and model compositions of the coexisting phases as additional constraints.

The models described in this study rely on a detailed characterization of the compositions and texture of the studied samples. Standardized X-ray maps are used to constrain the average composition of each growth zone and to calculate the phase proportions.

The GrtMod program was successfully used to model garnet growth conditions of a poly-metamorphic micaschist from the Sesia Zone (Western Alps). Garnet core records Permian HT/LP metamorphic conditions, whereas rims are formed at HP and MP during the Alpine continental subduction.

Acknowledgements

The authors acknowledge fruitful discussions with Rob Berman, Christian de Capitani, A. Pourteau, S. Centrella and V. Laurent on ideas regarding thermodynamic modeling and/or GrtMod. We also thank E. Duesterhoeft for adding new features to the program Theriak_D. We acknowledge H. Vrijmoed and an anonymous reviewer for their helpful comments and suggestions. This work was supported by the Swiss National Science Foundation (project 200020-146175).

Appendix 1 Automated strategy description

A1.1 Optimization1 (P–T)

Optimization1 is carried out within a P–T window defined between Tmin, Tmax, Pmin and Pmax (Fig. 4a, commands TMIN, TMAX, PMIN and PMAX see Appendix 2). The values of L0 and C0 are set at 1e19 outside this P–T window or if garnet is not stable. The minimum fraction of garnet to be stable during this stage is fixed by the estimated proportion in the present-day sample. During optimization1, garnet resorption is not allowed to change, and the volume fractions vi;j=1:i−1 of garnet stabilized during the previous stages are either entirely fractionated or partially fractionated from the bulk-rock composition. For advanced stages, the complete fractionation may generate extreme LEB compositions from which garnet become unstable in the Gibbs free energy minimization (Konrad-Schmolke et al., 2008). Four initial guesses are defined, as illustrated in Fig. 4 ( G1,1i, G1,2i, G1,3i, G1,4i). Additional initial guesses (G1,r>4i) can be easily defined. The number of initial guesses defines how many minimizations are done.

For stage S1 there are only two variables P1 and T1, hence optimization2 is skipped. At the end of each minimization a new solution is defined if (1) C0 is lower than STOL (see Appendix 1) and (2) if no previous solution exists with similar P and T (Fig. 4b). The P–T couple is not stored as a new solution if pressure and temperature differences with existing solutions are within TDI1 and PDI1, respectively (see Appendix 2). In this case the program considers that both minimizations converged to the same minimum, and only the first is stored as a solution. By contrast, for stage Si > 1, all solutions are stored, and the P–T couple with the smaller value of C0 is selected to be used as starting guess during optimization2 (for example S1,2i becoming G2,1i in Fig. 4b).

A1.2 Optimization2 (P–T–X)

Optimization2 is carried out for stage Si > 1 and the following ones. It involves i−1 additional compositional variables vi;=1:i−1 corresponding to the volume fractions of previous garnet growth zones. Pressure and temperature conditions of the best solution (S1,besti) from optimization1 are selected as initial guess (G2,besti). By contrast to optimization1, compositional variables allowing garnet resorption are introduced. The first initial guess is the exact solution of optimization1, without any resorption or with a fixed amount of resorption (Initial vi,j → max vj)). The first minimization allows testing if resorption can help to get a smaller value of C0 and therefore improve the quality of the solution. The second guess assumes very major resorption of previous garnet growth zones (Initial vi;jmin(vj)) and the third one an intermediate resorption (Initialvi,j=min(vj)+max(vj)min(vj)2). At the end of each minimization a new solution is defined if C0 is lower than STOL. The refinement phase is applied to all solutions.

A1.3 Auto-refinement phase

The auto-refinement phase aims to explore the P–T local variability of the cost function in order to provide a relative uncertainty and to investigate the sensitivity of the model compositions. The C0 value of the cost function is used because it intuitively represents the deviation between model and measured compositions. New C0 values are iteratively computed around the solution across height directions (D1, D2, …, D8 in Fig. 5). The P–T increments dT and dP are set using values defined in TDI1 and PDI1 (Appendix 1).

For a given solution S2,si with a value Co2,si and a tolerance TC0 (defined in RESC, see Appendix 1), across the direction d, the next P–T point n + 1 is calculated while the following criterion is met  
Con+1d<Co2,si+TCo.
(A1:1)

The same procedure is repeated in all directions in order to derive uncertainty bars (Fig. 5). The uncertainty on the volume fraction of garnet stable at Pi,Ti is estimated as the standard deviation of the volume fractions estimated at each point across all directions.

A1.4 Go fast mode

The go fast mode allows beginning optimization2 from a different starting point in order to check for alternative solutions. User defines the initial P–T couple and the program skips optimization1. As described in Appendix A1.2, three initial guesses with different compositional variables vi;j=1:i−1 (no resorption, strong resorption and moderate resorption) are defined for optimization2.

Appendix 2 GRTMOD main commands