An inverse modeling approach to obtain P–T conditions of metamorphic stages involving garnet growth and resorption

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. Key-words: thermodynamic modeling; garnet; resorption; X-ray mapping; XMapTOOLS.


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 (C BR ), 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 C BR , 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., 1984Spear et al., , 1991aKonrad-Schmolke et al., 2005;Gaidies et al., 2006Gaidies et al., , 2008bCheng & 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 C BR 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 C BR . -The second approach consists of successive forward models for which C BR 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 C BR , 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-Tstages 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).

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 S i (i > 1), the volume fractions of all previous garnet growth zones that are fractionated from the C BR . 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.

Growth stages of garnet and corresponding variables
As discussed in the introduction, in low-to moderategrade 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).
In the present study, the growth history of garnet is divided into discrete events, defined as "growth stages" (S i , i = 1,...,n, see Table 1). An individual growth stage is a short event occurring at given P i and T i during which a garnet volume fraction v i grows with a homogeneous composition C i grt . For any growth stage S i , the variables to be optimized are its specific P and T (P i and T i ) as well as the volume fraction of previous garnet to be fractionated from the C BR (v i,j = 1 ,..., v i,j = iÀ1 ). The number of variables at stage i thus is: For the first stage, the C BR of the sample is assumed to be equal to the composition of the reactive part of the system; hence C BR serves as input for the forward models. For subsequent stages (i > 1), the local effective bulk composition, C i LEB , is calculated as follows: where C i grt is the measured composition (in oxide wt%) of the garnet growth zone i; r i grt is its model density, and r iÀ1 rock and r iÀ1 mtx are the average rock density and matrix density (i.e., all phases except garnet) from the preceding stage; v i,j is the volume fraction of garnet crystallized during stage j that is fractionated from the C BR 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: The resorption of garnet j at stage i (i > j) is expressed in vol% of garnet in the bulk rock and is given by: The total resorption of garnet r grt is defined as In order to compute the local effective bulk composition from garnet volume fractionsto fractionate the previous garnet from the C BRthe rock density of the previous stage is required (Eq. (2)). The density of the rock at stage iÀ1 (r iÀ1 rock ) is calculated assuming zero porosity:  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).

Objective function
The objective function is composed of a loss function generating the number L 0 and a cost function generating the number C 0 . 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: where f measured k and f modeled k are the proportions of endmember k, as measured or modeled, and v k is the weighting factor of the corresponding element. The end-member proportions of grossular (f Grs ), pyrope (f Prp ), almandine (f Alm ), and spessartine (f Sps ) are calculated from Ca, Mg, Fe, and Mn abundances, expressed in number of atoms per formula unit (apfu). However, the proportions f Grs , f Prp , f Alm and f Sps 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 s 2 k (s k : standard deviation), hence the weighting factor is defined as In mapping conditions, the precision (1s) of the electron microprobe measurement of one pixel composition is estimated using a Poisson law (e.g., Lanari et al., 2014) 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 f k of the same X-ray flux f, the variance is close to  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 materialsuch as a set of pixels of X-ray mapsthe 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 n k of the corresponding element using the relationship Consequently, the loss function used in this study is The number L 0 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 C 0 , which is intuitively representing the quality of the solution. The cost function is defined as: C 0 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 C 0 < 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 C 0 (see application example below).

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.

Toward an automated optimization strategy
The optimization strategy depends on the number of variables being optimized. For instance, for the first stage S 1 , only two variables (P 1 , T 1 ) are optimized, and the procedure is straightforward. For all subsequent stages (S i > 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 (v i,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 (C 0 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.

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., 2012Lanari et al., , 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.

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 SiO 2 -Al 2 O 3 -FeO-MnO-MgO-CaO-Na 2 O-K 2 O-H 2 O using the same thermodynamic dataset (TC321p2.txt) and C BR 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 mm; Fig Grt 1 (0.37 vol%, Alm 36 Prp 2 Grs 26 Sps 36 ); Grt 2 (0.12 vol%, Alm 57 Prp 3 Grs 33 Sps 7 ); Grt 3 (4.2 vol%, Alm 67 Prp 7 Grs 25 Sps 1 ); and Grt 4 (2.31 vol%, Alm 64 Prp 9 Grs 26 Sps 1 ). 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 Fig. 4. Sketches illustrating the automated procedure. (a) Pressuretemperature diagram within T min -T max , P min -P max . The G i 1;r 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 C 0 of S i 1;2 is smaller than C 0 of S i 1;1 and is selected as the best P-T couple obtained during optimization1.  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 (C 0 between 0.007 and 0.025) resulting in an excellent match of the modeled compositions.

Sample description and compositional mapping
The studied sample FG12-157 is an eclogitic garnetbearing 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 mm 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 mm). 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. (2013Lanari et al. ( , 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 (SiO 2 ), anorthite (Al 2 O 3 , CaO), almandine (FeO), spinel (MgO), orthoclase (K 2 O), albite (Na 2 O), ilmenite (TiO 2 ), and tephroite (MnO). Analytical conditions for Xray 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 f Prp , f Grs and f Alm is about 0.03. This value is used to define the value of STOL.

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 (Alm 68-70 Prp 24-26 Grs 3-5 Sps 1-3 ,  (Fig. 1). Two distinct overgrowths surround the apparently porphyroclastic cores: a first rim (Alm 63-65 Prp 18-20 Grs 15-17 Sps 1 ) and a second rim (Alm 57-61 Prp 14-16 Grs 24-28 Sps 1 ) (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 (Grt 1 ), stage 2 for the first rim (Grt 2 ), and stage 3 for the second rim (Grt 3 ).

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  (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 Grt 1 , 3 vol% for Grt 2 , and 4 vol% for Grt 3 . 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.
The first minimization (G 1 1;1 from 600°C-8 kbar) converges to a minimum at 851°C and 6.03 kbar for a C 0 value of 0.021. A solution S 1 1;1 is retained for the following because for this first minimization C 0 < 0.03.
Model f Alm and f Grs 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 (G 1 1;2 from 600°C-16 kbar) converges to a different minimum at 674°C and 17.49 kbar with a C 0 value of 0.107. As C 0 is much higher, model f Alm and f Grs 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 C 0 > 0.03. The third minimization (G 1 1;3 from 800°C-8 kbar) converges to a minimum at 851°C and 6.02 kbar for a C 0 value of 0.021. This result is fairly similar to solution S 1 1;1 , indeed minimizations 1 and 3 converge to the same local minimum. The fourth minimization (G 1 1;4 from 800°C-16 kbar) converges to a minimum at 899°C and 6.14 kbar for a C 0 value of 0.015. A solution S 1 1;2 is obtained (C 0 < 0.03). 7.14 vol% of garnet is predicted to be stable at that stage. The second solution (S 1 1;2 ) has a smaller value of C 0 (Co 1 1;2 < Co 1 1;1 ) and is selected as the best solution for stage 1: S 1 best ¼ S 1 1;2 (see Fig. 7; Tables 3 and 4).

Stage 2automated procedure
For stage 2, the optimization is divided into two steps: optimization1 and optimization2. During optimization1, a fixed amount of garnet Grt 1 (7.14 vol%), corresponding to the amount predicted during stage 1, is initially fractionated from the C BR . The P-T conditions of the four starting guesses G 2 1;1 , G 2 1;2 , G 2 1;3 , G 2 1;4 are the same as for stage 1 (see above). The first minimization (G 2 1;1 from 600°C-8 kbar) converges to a minimum at 719°C and 9.66 kbar with C 0 = 1.74 Â 10 À4 . A temporary solution S 2 1;1 is stored; 5.95 vol% of garnet is predicted to be stable at that stage. The second minimization (G 2 1;2 from 600°C-16 kbar) converges to a minimum at 647°C and 16.03 kbar with C 0 = 4.15 Â 10 À3 . A temporary solution S 2 1;2 is stored; 9.02 vol% of garnet is predicted to be stable at that stage. The third minimization (G 2 1;3 from 800°C-8 kbar) converges to a solution S 2 1;3 at 719°C and 9.66 kbar with C 0 = 1.25 Â 10 À4 . This solution is close to S 2 1;2 with a slightly smaller residual. The fourth minimization (G 2 1;4 from 800°C-16 kbar) converges to a solution S 2 1;4 similar to S 2 1;2 . Optimization1 of stage 2 shows that for the same C BR (i.e., without resorption) Grt 2 is predicted stable at 719°C and 9.66 kbar and 647°C and 16.03 kbar. The automated algorithm selects S 2 1;3 as the best solution based on the C 0 values (S 2 1;best ¼ S 2 1;3 ). The P-T conditions of initial guesses of optimization2 are fixed at 719°C and 9.66 kbar. A new variable v 2,1 corresponding to the quantity of garnet Grt 1 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 (v 2,1 = 7.14 vol%), strong resorption (v 2,1 = 4 vol%) and moderate resorption (v 2,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 (G 2 2;1 of 7.14 vol%) converges to a minimum at 719°C and 9.65 kbar and a final v 2,1 of 7.13 vol%, with a C 0 value of 6.55 Â 10 À5 and 5.95 vol% of Grt 2 . The second minimization (G 2 2;2 of 4 vol%) converges to a minimum at 694°C and 9.49 kbar and a final v 2,1 of 4.35 vol% with a C 0 value of 4.44 Â 10 À5 and 9.44 vol% of Grt 2 . The third minimization (G 2 2;3 of 5.57 vol%) converges to a minimum at 707°C and 9.57 kbar and a final v 2,1 of 5.86 vol% with a C 0 value of 2.46 Â 10 À5 and 7.53 vol% of Grt 2 . These results suggest that Grt 2 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 Grt 2 (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.

Stage 2go 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 S 2 1;2 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 v 2,1 of 5.92 vol% with C 0 = 8.12 Â 10 À6 and 10.40 vol% of Grt 2. 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 (DT = 0.37°C, DP = 0.09 kbar). In this case, resorption of 1.21 vol% of Grt 1 improves the quality of the model and is selected as the most likely solution for stage 2.

Stage 3automated procedure
For stage 3, the optimization is again divided into two steps: optimization1 and optimization2. During optimi-zation1, a fixed quantity of garnet Grt 1 and Grt 2 (5.92 vol% and 6.00 vol%) is fractionated from the C BR . For Grt 2 optimization1 is carried out assuming 4.40 vol% resorption, the minimum required to stabilize garnet with a composition similar to Grt 3 . The first minimization (G 3 1;1 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 Grt 3 is fixed at 4 vol% observed in the sample). The second minimization (G 3 1;2 from 600°C-16 kbar) converges to a minimum at 637°C and 13.33 kbar for C 0 = 1.71 Â 10 À4 . A temporary solution S 3 1;1 is stored; 6.17 vol% of garnet is predicted stable at that stage. The third minimization (G 3 1;3 from 800°C-8 kbar) does not converge because it predicts only 2.2 vol% of garnet, far below the minimum amount for  Grt 3 . The fourth minimization converges to the same minimum as S 3 1;1 with a slightly smaller C 0 value of 1.62 Â 10 À4 . A temporary solution S 3 1;2 is stored. The algorithm selects S 3 1;2 as the best solution based on the C 0 values (S 3 1;best ¼ S 3 1;2 ). The P-T conditions of initial guesses of optimization2 are fixed at 637°C and 13.33 kbar. Two new variables v 3,1 and v 3,2 are to be optimized; they correspond to the quantities of garnet Grt 1 and Grt 2 fractionated from the bulk composition during stage 3. Three starting guesses are defined assuming either no more resorption than given by the input value (v 3,1 = 5.92 vol% and v 3,2 = 6.00 vol%), strong resorption (v 3,1 = 4.00 vol% and v 3,2 = 3.00 vol%), and moderate resorption (v 3,1 = 4.96 vol% and v 3,2 = 4.50 vol%). In this case initial v 3,2 is set at 6.00 vol% whereas the volume fraction of Grt 2 at stage 2 is 10.40 vol%. As only a small volume fraction of Grt 2 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 (S 3 2;1 ) at 637°C and 13.33 kbar and v 3,1 and v 3,2 of 5.92 and 6.00 vol% with C 0 = 9.72 Â 10 À5 and 6.16 vol% Grt 3 . The second minimization converges to a minimum (S 3 2;2 ) at 632°C and 11.01 kbar and v 3,1 and v 3,2 of 4.86 and 3.40 vol % with C 0 = 2.68 Â 10 À6 and 8.42 vol% of Grt 3 . The third minimization converges to a minimum (S 3 2;3 ) at 637°C and 13.33 kbar and v 3,1 and v 3,2 of 5.16 and 4.62 vol% with C 0 = 5.89 Â 10 À6 and 8.27 vol% Grt 3 . Two solutions are found at slightly different pressures (13.37 kbar and 11.01 kbar). The second solution S 3 2;2 is selected here because it is considered more likely based on the C 0 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.

Intergranular diffusion and global equilibrium within domains
In this study garnet growth is modeled based on C BR 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).

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 C 0 and L 0 values for a given effective bulk composition is provided in GRTMOD.
As an example, the C 0 map in the P-T range 500-900°C and 5-20 kbar was computed for stage 1 using the C BR (Fig. 8a,b). This map exhibits the shape of the cost function C for stage 1 (see Sect. 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, C 0 = 0.107) and the second region at high T (HT; 850-900°C, ∼6 kbar, C 0 < 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.

Automated strategy [1]: limitation of multiple minima and solution finding
Similarly, the C 0 map in the range 500-900°C and 5-20 kbar was computed for stage 2, assuming no resorption of Grt 1 (Fig. 8c,d). The effective C BR is computed by subtracting from C BR the amount of Grt 1 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 C 0 of 1.25 Â 10 À4 ), and the other at HP (647°C and 16 kbar with C 0 = 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 C 0 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 (C 0 = 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 C 0 between both solutions, and Grt 2 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 Grt 2 suggest that it coexisted with Si-rich phengite (Si 4þ = 3.33 ± 0.02 apfu, X Mg = 0.76 ± 0.02) ( Table 5). The K-rich white-mica composition predicted by the model at 647°C/16 kbar is Si 4þ = 3.34 and X Mg = 0.79, whereas at 720°C/6.6 kbar it is Si 4þ = 3.17 and X Mg = 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 Grt 2 growth during a HP stage.

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. The C 0 maps in the P-Trange 550-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.

P-T stages recorded in garnet from a polymetamorphic micaschist
The model predicts that garnet core (Grt 1 ) 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 (Grt 2 ) 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 Grt 2 to grow at the expense of Grt 1 . The shape of Grt 1 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 Grt 1 is very different from any observed in comparable rocks containing typical Alpine prograde garnet. The P-T conditions obtained for Grt 1 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 Grt 1 and the precipitation of Grt 2 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 (Grt 3 ) 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.

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 Grt 2 was resorbed. Such partial resorption has a strong effect on the effective C BR 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 Grt 1 was initially fractionated from the C BR 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. Garnet resoption (in % of garnet core) 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.) -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 (P 1 -T 1 ) are found, then a second point is analyzed and finds P 2 -T 2 , 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.  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 Grt 2 , Grt 3 and Grt 4 ). 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 Grt 2 and Grt 3 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-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 Grt 2 , Grt 3 and Grt 4 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-1 and MP-2) fail to reproduce the observed zoning profile. Furthermore, the fraction of Grt 2 is always overestimated at the first stage of growth (∼16 kbar). For Grt 3 , 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 (Grt 3 and Grt 4 in Fig. 10c). In contrast, model MP-2 does match the composition of Grt 3 but the corresponding volume fraction is seriously underestimated, whereas Grt 2 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.

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. initial guesses (G i 1;r>4 ) can be easily defined. The number of initial guesses defines how many minimizations are done.
For stage S 1 there are only two variables P 1 and T 1 , hence optimization2 is skipped. At the end of each minimization a new solution is defined if (1) C 0 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 S i > 1, all solutions are stored, and the P-T couple with the smaller value of C 0 is selected to be used as starting guess during optimization2 (for example S i 1;2 becoming G i 2;1 in Fig. 4b).

A1.2 Optimization2 (P-T-X)
Optimization2 is carried out for stage S i > 1 and the following ones. It involves iÀ1 additional compositional variables v i;j¼1:iÀ1 corresponding to the volume fractions of previous garnet growth zones. Pressure and temperature conditions of the best solution (S i 1;best ) from optimization1 are selected as initial guess (G i 2;best ). 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 v i;j → maxðv j Þ).
The first minimization allows testing if resorption can help to get a smaller value of C 0 and therefore improve the quality of the solution. The second guess assumes very major resorption of previous garnet growth zones (Initial v i;j → minðv j Þ) and the third one an intermediate resorption (Initial v i;j ¼ minðv j Þ þ maxðv j ÞÀminðv j Þ 2 ). At the end of each minimization a new solution is defined if C 0 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 C 0 value of the cost function is used because it intuitively represents the deviation between model and measured compositions. New C 0 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 S i 2;s with a value Co i 2;s and a tolerance TC 0 (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 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 P i ,T i 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 v i;j¼1:iÀ1 (no resorption, strong resorption and moderate resorption) are defined for optimization2. T steps first level TDI2 T steps second level PDI1 P steps first level PDI2 P steps second level