An Internally-Consistent Database for Oxygen Isotope Fractionation Between Minerals

The knowledge of the fractionation behaviour between phases in isotopic equilibrium and its evolution with temperature is fundamental to assist the petrological interpretation of measured oxygen isotope compositions. We report a comprehensive and updated internally consistent database for oxygen isotope fractionation. Internal consistency is of particular importance for applications of oxygen isotope fractionation that consider mineral assemblages rather than individual mineral couples. The database DBOXYGEN is constructed from a large dataset of published experimental, semi-empirical and natural data, which were weighted according to type. It includes fractionation factors for 153 major and accessory mineral phases and a pure H2O fluid phase in the temperature range of 0–900 C, with application recommended for temperatures of 200–900 C. Multiple primary data for each mineral couple were discretized and fitted to a model fractionation function. Consistency between the models for each mineral couple was achieved by simultaneous least square regression. Minimum absolute uncertainties based on the spread of the available data were calculated for each fractionation factor using a Monte Carlo sampling technique. The accuracy of the derived database is assessed by comparisons with previous oxygen isotope fractionation calculations based on selected mineral/mineral couples. This database provides an updated internally consistent tool for geochemical modelling based on a large set of primary data and including uncertainties. For an effective use of the database for thermometry and uncertainty calculation we provide a MATLABVC -based software THERMOOX. The new database supports isotopic modelling in a thermodynamic framework to predict the evolution of dO in minerals during metamorphism.


INTRODUCTION
Stable isotopes are important tools for a wide range of applications in Earth Sciences as the isotopic composition of minerals can record their physical and chemical conditions of equilibration. Oxygen isotope fractionation between two cogenetic minerals is, for example, temperature-dependent and has been intensively used as mineral thermometer (e.g. Bottinga & Javoy, 1973;Clayton, 1981;Zheng, 1993a). Additionally, the oxygen isotope composition of co-existing phases is a prime tool for reconstructing fluid-rock interaction and evaluating mineral equilibration (e.g. Valley, 1986Valley, , 2001Eiler et al., 1992Eiler et al., , 1993Kohn, 1993;Valley & Graham, 1996;Baumgartner & Valley, 2001;Page et al., 2014;Rubatto & Angiboust, 2015). Such applications require the knowledge of the fractionation behaviour between two mineral phases as well as the possible evolution with temperature. Over the past decades, efforts have been directed towards the determination of fractionation factors between minerals or minerals and fluids. Pioneer works were mostly based on experimental data on isotope exchange between minerals and water Clayton et al., 1972;Matsuhisa et al., 1979;Matthews et al., 1983a), whereas mineral/mineral couples were often derived by combining these calibration data.
Phase equilibrium modelling of metamorphic rocks has evolved to a point where mineral assemblages are being considered, rather than individual mineral couples (see Lanari & Duesterhoeft, 2019 for a review). A step forward can also be anticipated for oxygen thermometry and isotopic tracing, provided that the fractionation factors between several minerals can be combined in a similar way to that done for the thermodynamic data of minerals (e.g. Berman, 1988;Holland & Powell, 1998) or aqueous species (Miron et al., 2016(Miron et al., , 2017. Ideally, the achievement of internal consistency can be defined as the optimal satisfaction of the simultaneous constraints imposed by all the data. The internal consistency of a database can be formulated as a list of increasingly stringent conditions (Engi, 1992), including: (1) compatibility with theoretical definitions and reference values; (2) representation of all of the listed data, in the sense of 'best compromise' with all of the available data considered simultaneously; and (3) compatibility with all underlying data within stated uncertainties.
For oxygen isotope fractionation, early works combining experimental data on different minerals resulted in self-consistent compilations (e.g. Bottinga & Javoy, 1975;Chiba et al., 1989;, each of them involving typically less than 10 minerals. The more recent compilation by Chacko et al. (2001) based on carbonate exchange experiments includes 15 minerals. A more considerable effort is needed for critical assessment of a larger pool of primary data and for their transformation into an internally consistent dataset. A second category of database for oxygen isotope fractionation factors is based on semi-empirical calculations such as the theory of bond strength models (Garlick, 1966;Schü tze, 1980). For instance, the work of Zheng (1991Zheng ( , 1992Zheng ( , 1993aZheng ( , 1993bZheng ( , 1993cZheng ( , 1996 resulted in a large database of semi-empirical data and offering partial compatibility within individual mineral groups . In that approach, the change of reference mineral from quartz (Zheng, 1991(Zheng, , 1992(Zheng, , 1993a(Zheng, , 1993b to calcite (Zheng, 1996 does not guarantee consistency between the different groups. Internal consistency is a fundamental criterion for advanced applications of oxygen isotope fractionation that consider mineral assemblages rather than individual mineral couples. However, internal consistency is difficult to achieve as there are several reasons, such as experimental errors, inadequate theory or assumptions, or poor measurements, that may cause minor or major discrepancies between different calibrations (May & Murray, 2001). In this regard, combination of natural, experimental, theoretical and semi-empirical data could be used to overcome the limits of each method. This paper presents a large and internally consistent database for oxygen isotope fractionation between minerals. The first section describes how existing oxygen isotope fractionation data were selected, weighted and processed. We offer a different approach than other databases that provide calculations based on single calibration chosen by the user (see AlphaDelta by G. Beaudoin & P. Therrien available at http://www2.ggl. ulaval.ca/cgi-bin/alphadelta/alphadelta_4alpha.cgi). Our approach considers simultaneously a large number of mineral couples and maximises the consistency of the oxygen isotope fractionation among all of them. It enables calculations to be done using phase couples for which no direct data is available, and it is directed to modelling rock systems where the 18 O/ 16 O composition of every phase is required to be in equilibrium with every other phase, in an internally consistent manner. The significance of the results and their uncertainties is discussed in detail and possible applications in metamorphic petrology are shown. The internally consistent database DBOXYGEN version 2.0.3 contains fractionation factors for 153 minerals and a pure water fluid phase based on a large compilation of literature data (Table 1). We also present the MATLABV C -based software THERMOOX that uses this database for thermometry and uncertainty calculation. The database, supporting information, software solutions and further updates are available for download at http://oxygen.petrochronol ogy.org/.

Methods for estimating oxygen isotope fractionation between phases
Controlled laboratory experiments using natural or synthetic phases are the most commonly used method for estimating oxygen isotope fractionation between mineral and fluid phases. The fractionation factors are obtained by interpolating the experimental data points along the temperature range [T min , T max ] of the experiments. However, the method of equilibrium exchange requires temperatures where isotopic exchange is sufficiently rapid to reach equilibrium in a reasonable amount of time, typically T > 500 C. The success of this type of experiments can be limited by small exchange amounts and failure to attain complete isotopic equilibrium (O'Neil, 1986;Chacko et al., 2001). Significant advances in experimental techniques are represented for instance by (1) the development of the three isotope exchange method to estimate equilibrium fractionation from partially exchanged samples (e.g. Matsuhisa et al., 1978;Matthews et al., 1983aMatthews et al., , 1983bMatthews et al., , 1983c and (2) the introduction of piston-cylinder high-pressure exchange experiments that allow mineral-water and mineralcalcite exchange to occur much faster than in the previously used conventional 1-2 kbar hydrothermal apparatus (e.g. Clayton et al., 1975;Matsuhisa et al., 1978;Matthews et al., 1983aMatthews et al., , 1983bMatthews et al., , 1983c. The compilation by Chacko et al. (2001) revisited in detail experiments published before 2001.
The third approach to obtain fractionation factors is the investigation of natural samples (e.g. Bottinga & Javoy, 1975;Valley, 2003;Lacroix & Vennemann, 2015). The analysed mineral phases are assumed to have coexisted in equilibrium at a temperature that is determined with an independent thermometer and the isotopic compositions to be preserved during cooling. In general, experiments and semi-empirical calculations have been favoured, largely due to the recognition of widespread diffusional resetting of isotope in natural samples, making many natural rocks inappropriate for retrieving equilibrium fractionation factors (e.g. Giletti, 1986;Eiler et al., 1992Eiler et al., , 1993Valley, 2001).
If the fractionation parameters between two couples i-j and i-k are known, parameters for the combined couple k-j can be obtained from the algebraic sum

Minerals involving solid solutions
It has been demonstrated that cation substitution can influence significantly the fractionation between minerals Zheng, 1993aZheng, , 1993bKohn & Valley, 1998b. Available experimental data for feldspars and muscovite  and thermodynamic calculations for plagioclase (Ganguly, 1982) show that solid solutions are well approximated by linear interpolation between endmembers. Although the effect of cation substitution on 18 O/ 16 O fractionation has been investigated via experiments only for few silicate groups, these relationships are commonly accepted and applied (e.g. Kohn, 1993;Kohn & Valley, 1998b). Hence the isotopic composition of a phase i involving a solid solution with k phase components (end-members) can be recalculated as a linear combination of the oxygen isotope compositions of each end-member d 18 O kÀi (Kohn, 1993) with X kÀi the molar fraction of end-member k in phase i, N kÀi and N i the number of moles of oxygen in the endmember k and in phase i. The isotope fractionation parameters between a solid solution i involving k phase components and a pure phase j are obtained from By combining Equation (4) with Equations (8-10) and assuming the validity of the approximation in Equation (5), the stable isotope fractionation between two phases i (with k phase components) and j (pure) as a function of T becomes: Equation (11) can be rearranged to obtain a set of linear equations valid for each phase component.

Effect of pressure and salinity in experimental settings
The effect of pressure on oxygen isotope fractionation between minerals has been demonstrated to be small (typically < 0Á2&), particularly at high temperature (Clayton et al., 1975;Polyakov & Kharlashina, 1994;Chacko et al., 2001). Our calculation according to the method of Polyakov & Kharlashina (1994) of pressureinduced deviations between 0Á1 and 1Á0 GPa at different temperatures for key mineral couples involving quartz confirm this ( Table 2). The deviations are higher for quartz/calcite and quartz/rutile than for the silicate/silicate couples, but they are far below typical uncertainties of the fractionation models, especially at T ! 350 C. Therefore, the effect of pressure is not considered for the calculation of our database. Truesdell (1974) investigated the effect of dissolved salts on the oxygen isotope activity ratio of water up to 275 C. The author concluded that most of the salt solutions show a complex behaviour of oxygen isotopes with increasing temperature and that mineral-water fractionations may be affected by shifts up to $ 3& in both the directions. However, this observation is in contrast with other studies investigating different chemical systems. Horita (1989) reviewed 'salt effect' coefficients measured by several authors and investigated the effect of high-CaCl 2 or -MgCl 2 brines on hydrogen and oxygen isotope fractionation factors. The author concluded that NaCl, as well as gypsum, has not effect, although most data were obtained at low temperature (i.e. < 40 C). According to the experiments of Kendall (1983), the presence of NaCl up to 4 molal concentration has no effect on oxygen isotope fractionation between calcite and water at 275 C. Zhang et al. (1989) studied experimentally oxygen isotope fractionation between quartz and water, and also concluded that the effect of salinity is not significant at temperature higher than 250 C, whereas it probably affects the rate of oxygen isotope exchange among oxygen-bearing phases. Zhang et al. (1994) observed the same effect for the quartzcassiterite-water and the wolframite-water systems. The experiments performed by  show instead that, at elevated temperature and pressure, mineral/mineral fractionation factors derived from separate mineral-'pure' water experiments are not fully compatible because of the different species dissolved in the water in each experiment.
In our database, water is always assumed to be a pure H 2 O phase, and experiments performed under different NaCl salinities are still considered. However, because previous studies provide reasons to be cautious about experiments that used brines, such experiments are weighted less than primary data using pure water (see below). Details on the use of brines in experiments are reported in Supplementary Data 1, available for downloading at http://www.petrology.oxfordjournals. org.

Selection of primary data and temperature range
In this study, a wide compilation of published fractionation factors calculated by the methods mentioned above has been produced. Experimental, semiempirical and natural data (primary data, Fig. 1) have been compiled from the literature for 400 mineral/aqueous fluid and mineral/mineral couples including most major and accessory phases. The complete list of phases is reported in Table 1, the list of all the phase couples in Supplementary Data 1 and 3. All the available fractionation functions between two pure phases or phase end-members were used for the calculation, with the exceptions of calibrations that have been corrected or superseded by later studies, in which case All the available oxygen isotope fractionation factors (experimental data, natural studies and semi-empirical calculations) for a phase couple (primary data) have been considered when not inconsistent with the other available data (see text for details). (b) Secondary data points have been generated for each primary dataset by discretizing the fractionation function along the temperature range [T min , T max ] with a temperature step T inc . This parameter is used as a weighting factor and reflects the degree of confidence of the corresponding primary data. A fixed uncertainty of 6 0Á3& was attributed to each secondary data point (not shown in the plot, see text for details). The red line represents the simple interpolation model (second order polynomial fit) to all the secondary data for each phase couple. The results (best fit) of this iterative interpolation are not internally consistent. They represent a first approximation used as input for the global optimization that delivers internal consistency. Table 2: Pressure-induced corrections to the fractionation values at P ¼ 1 GPa relative to ambient pressure for chosen mineral couples. Calculation by using the method of Polyakov & Kharlashina (1994). Mineral abbreviations from Whitney & Evans (2010) only the revised values are reported (i.e. Northrop & Clayton, 1966;Clayton et al., 1972revised by Friedman & O'Neil, 1977 partially revised by Zheng, 1996;Zheng & Bö ttcher, 2016 partially revised by Rubatto et al., 2014 and and of calibrations that have been demonstrated to be unreliable (e.g. Shiro & Sakai, 1972;Krylov et al., 2002;Zhou & Zheng, 2003, 2005. Calibrations that have been evaluated but not used for the final calculation are reported in Table 1 in italics and details are given in Supplementary Data 1. Fractionation factors for solid solutions calculated without distinguishing the contribution of each phase component (e.g. Javoy et al., 1970;Bottinga & Javoy, 1975;Kohn & Valley, 1998b) were mostly excluded. In such cases, the fractionation factors are valid only for specific mineral compositions and they can only be used to test the final database. As discussed above, the fractionation between two phases is generally reported in the literature as a function of T with three parameters A, B and C defined for a temperature range [T min , T max ] in which the experiments or the calculation were performed. For the cases where only experimental points were available, A, B and C were calculated by second order polynomial fitting (e.g. Becker and Clayton, 1976;. Bond-strength models (e.g. Richter & Hoernes, 1988;Zheng, 1991Zheng, , 1993aZheng, , 1993b present equations for calculating fractionation factors between phases over a wide range of temperature (e.g. 0-1200 C), whereas experimental and natural data never exceed the stability region of a phase and are usually even narrower. The temperature range selected for the presented database is between 0 C and 900 C (see discussion below), corresponding to the range where the most experiments have been performed in order to decrease the otherwise overemphasized contribution of the semiempirical calibrations. None of the primary data from experimental results or natural samples were extrapolated outside the prescribed temperature range, as the parameters are not always valid over a broad temperature range (Clayton, 1981). All the primary data published to be valid outside the chosen temperature range of 0-900 C have been excluded from the calculation or restricted to such range (see Supplementary Data 1 for details).
Theoretically, at infinite T the fractionation between any set of two phases approaches 0& and, therefore, many authors questioned the validity of the parameter C in Equation (4) (Sharp, 2017). Whereas the database was optimized assuming C ¼ 0, no primary data were modified after the published (or revised) functions, because this option was considered as an excessive interference with the primary data.

Generation of secondary data and weighting procedure
The experimental data points used to derive the primary data and their analytical uncertainties are not systematically available. This prevented the original experimental data points (and their uncertainties) to be used for fitting. Instead, secondary data points were generated for each couple by discretizing the fractionation function along the temperature range [T min , T max ] with a temperature step T inc (Fig. 1). Over a given temperature range, the number of secondary Td 18 O points increases with decreasing T inc and thus, a function discretized by using a smaller value of T inc will have more weight during the optimization (see below) than a function discretized with a larger value of T inc . This parameter is used as a weighting factor applied to the secondary data reflecting the degree of confidence of the corresponding primary data. If not specified differently in Supplementary Data 1, a T inc value of 5 C was used for well-constrained experimental and natural data, 10 C for synthesis experimental data or experiments with salt added, and 20 C for semi-empirical data.
In addition, a fixed uncertainty of 0Á3& (2r) was attributed to each secondary data point, based on the average values proposed in the literature Savin & Epstein, 1970;Cole, 1985;Zheng, 1991Zheng, , 1993bVitali et al., 2000). The uncertainty on the temperature measurement during an experiment is rarely reported in literature and it does not apply to semi-empirical calculations; therefore, uncertainty in T was not taken into account.

Fitting procedure and optimization
Notwithstanding that fractionation factors for mineral couples may be most precise for single, specific mineral compositions, the presence of inconsistency when comparing fractionation factors from different compilations remains an issue. An internally consistent database aims to improve (1) the accuracy by considering numerous independent data and different data types, and (2) the versatility because it can be applied to mineral assemblages without carrying inconsistencies between mineral pairs. Data processing is required to achieve internal consistency (as defined in the introduction) and it can become prohibitive as large multicomponent systems are involved. Manual evaluations are accordingly confined to small chemical systems or can only be performed infrequently (May and Murray, 2001). In addition, the development of a computational procedure allows the database to be easily updated when new primary data become available.
The internally consistent dataset of fractionation parameters A, B and C (Eq. 4) and their uncertainties were obtained from the secondary data following a three-step procedure: (1) a first iterative approximation of the fractionation parameters; (2) a least-square global optimization; and (3) the calculation of uncertainties by mapping the uncertainty envelopes using a Monte-Carlo (MC) technique. A series of MATLAB V C -based programs were used; short descriptions of each part are given in the following.
Step 1: Iterative approximation of fractionation factors The first step consists of iteratively applying for each mineral/mineral or mineral/water couple a simple interpolation model to the available secondary data as defined above (Fig. 1). The detailed procedure is illustrated in Fig. SD4-1 (1) it can be considered a pure phase; (2) it is one of the most common rock-forming minerals, with a wide range of P-T stability; (3) it offers the largest collection of data, derived from different theoretical and experimental methods; and (4) it has been widely used in literature (i.e. proposed fractionation factors for different couples mineral/quartz). It is important to mention that geothermometry or fractionation modelling can be performed in quartz-absent systems using the internally consistent database (after step 2, see below), because every combination of fractionation parameters between all the phases available in the database can be recalculated (see Eq. 6).
Step 2: Global optimization using least squares Once an approximation of a first solution has been made, a selection criterion or objective function that can be minimized is defined in order to derive a unique and optimal solution (e.g. Berman et al., 1986). The objective function was defined as the sum of the square of the differences between the d 18 O value Z l and the modelled d 18 O value Y l of the secondary data point l.
The minimization of this objective function is a nonlinear optimization problem, which can be treated with the Nelder-Mead algorithm (Nelder & Mead, 1965). The mineral/mineral or mineral/water couples for which a larger number of secondary data points are available are predominant in constraining the final dataset. The weighing attributed to the different data types (see above) in the generation of secondary data steers the final model to be more heavily determined by experimental and natural data.
Step 3: Determination of minimum uncertainties Uncertainties on oxygen isotope fractionation parameters serve the important purpose of indicating relative errors, i.e. to know which parameters are better constrained and consistent with available data, and which are less constrained. Unfortunately, most of the studies reporting primary data do not report analytical uncertainties. This lack of constraints prevents the uncertainties to be assessed with the highest rigor. Instead, for each fractionation couple, a 2r uncertainty envelope was obtained containing at least 95% of the secondary data points used for fitting.
The regression models are polynomial in one variable (T), and their coefficients A, B and C are strongly correlated. Because of this correlation, the coefficients cannot be considered to have independent sources of error. As there exists no standard analytical formula to compute or propagate the uncertainty of polynomial models (see the Guide to the Expression of Uncertainty in Measurement, ed. Bureau International des Poids et Mesures, 2008), it can only be estimated via numerical simulation or linearization approximation. To define the 95% uncertainty envelope, a numerical procedure involving MC analysis was used to approximate the uncertainties on each parameter A, B and C. The procedure is illustrated in SD4-2 in Supplementary Data 4. For were generated using an uniform distribution over a range of 6 100% (or 6 1& when the value was lower than 0Á001) for r(A REF-i ) and r(B REF-i ) and over the range of 6 1& for the parameter C ¼ 0. This interval for C was chosen to contain the average (-0Á95&) and the median (0Á14&) of the C values of all the considered experimental and natural primary data. The uncertainties on the parameters of the couples that do not directly involve the reference mineral (REF) were obtained using the following relationships: Any set of uncertainties [r(A i-j ), r(B i-j ), r(C i-j )] calculated for the phase couple (i-j) bounds a region with upper and lower limit of 1000lna sðu;lÞ defined as The deviation of each secondary data point (coordinates: T p , 1000lna p ) from the optimized model (1000lna m ¼ f ðT Þ in T ¼ T p ), noted as DPM (Distance-Point Model), was calculated and compared with the vertical distance (i.e. at fixed T) between the envelope 1000lna s and the model 1000lna m , noted as DEM (Distance-Envelope Model) (Fig. 2) , the one minimizing DEM over the whole T range was selected as final set of uncertainties (Fig. 2). The absolute values of the uncertainties (at 2r) are reported in Table 1 and the values of the correspondent residual (Eq. 12) in Supplementary Data 2.
It is important to note that the calculated values of r represent a minimum absolute uncertainty that is based exclusively on the spread of the available data. The position of each envelope depends on the number and on the spread of the secondary data generated. In cases where fractionation factors are constrained by a single set of secondary data, the calculated uncertainties are expected to be small due to the absence of spread in the data and they might increase when new secondary data are added. Thus, it is important to stress that small uncertainties may be indicative of a lack of constraints rather than of a good agreement among them. Calculated models for couples based on several sets of primary data (see Table 1) are expected to be more accurate, but they might show larger uncertainties due to the distribution of the data. The user of the database is invited to check Table 1 and Supplementary Data 1 and 3 in order to evaluate on how many sets of primary and thus secondary data the reported uncertainties are based.

MODEL RESULTS COMPARED TO PRIMARY DATA
The number of mineral phases and end-members of solid solutions considered in this study (Table 1) is controlled by the availability of oxygen isotope fractionation data. The database (Supplementary Data 2) includes a pure water phase, 25 groups of silicates and 13 single silicate minerals, for a total of 100 pure silicate phases and solid solution end-members, as well as 15 carbonates, 17 oxides, 12 hydroxides, 3 phosphates, 1 salt, and 4 sulphates. Table 1 presents the model results for the fractionation parameters A and B (C ¼ 0) (Eq. 4) between quartz and each phase after the optimization for internal consistency, and the calculated minimum absolute uncertainties at the 2r level (see above). Secondary data and models for key or exemplar phases are shown in Figs 4, 5 and 6. The secondary data and models for all the considered phase couples are reported in Supplementary Data 3.
The calculated fractionation parameters for silicates and carbonates are, in most cases, supported by a diversity of data from different techniques and thus considered to be more robust. Of the major rock-forming minerals, quartz, garnet and pyroxene are based on solid sets of primary data and are inferred to be well constrained. For solid solutions such as amphibole, micas, chlorite, and for minerals of complex chemical composition (e.g. epidote), the primary data are scarce or inconsistent and thus their fractionation factors may lack accuracy or bear a large uncertainty. Mineral groups such as spinel, ilmenite, humite, feldspathoid, pyroxenoid, phenacite, and most of the oxides and hydroxides are uniquely based on semi-empirical datasets by Zheng (1991Zheng ( , 1993aZheng ( , 1993bZheng ( , 1996 because of the lack of experimental and natural data. Experimental or natural data at T < 200 C are available for fractionations between quartz/water and involving carbonates, chlorite, serpentine, zeolites, few oxides and sulphates. Despite their paucity and possible spread, those data are of key importance beside semiempirical data for constraining the shape of the fractionation functions, and they improve the quality of the fit also at higher T (see next section).
Further details, addressing in particular the inconsistencies among literature data and the behaviour of the model for those cases are discussed below. Although the parameters for each mineral were derived simultaneously so that every dataset contributed to the refinement of the fractionation properties, the following description is organised by mineral groups. Results are discussed between 200 and 900 C, with emphasis on the stability range for each phase (Figs 3,4,5,6). Results at T < 200 C are discussed only for quartz/water and calcite/water pairs.

Anhydrous silicates Quartz
Quartz was chosen as reference mineral, i.e. the parameters A and B (C ¼ 0) in Table 1 describe the fractionation between quartz and a phase. The fractionation between quartz and water is the best documented in the literature, with 16 datasets (Table 1) derived with different methods, of which 14 were selected for the optimization. Their good consistency is demonstrated by the low uncertainties on the fit reported in Equation (17) (at 2r) and plotted in Fig. 3 The internally consistent optimized model (opt, Eq. 18, uncertainties at 2r) agrees with the fit of the quartz/water secondary data (fit, Eq. 17) within 0Á4&, proving good consistency not only among data for quartz/water, but also among all the data involving quartz (135 couples) and water (147 couples). The uncertainties on the fractionation parameters in Equation (18) incorporate not only the spread of secondary data for this couple, but also for all the other secondary data involving either quartz or water. Sharp et al. (2016) proposed a calibration of the triple oxygen isotope fractionation for quartz/water based on high-T exchange experiments and low-T empirical estimates. All the data considered by the authors are also included in our dataset independently. Our optimized model agrees within 0Á5& with the calibration proposed by Sharp et al. (2016) at T ! 200 C, and it is 1& lower at 60 C.
Theoretically, because mineral polymorphs have a different crystal structure, they are expected to behave differently with respect to oxygen stable isotope fractionation (Shiro & Sakai, 1972;Kawabe, 1978;Smyth, 1989;Zheng, 1993c). Nevertheless, the fractionation functions for a-quartz/water and b-quartz/water proposed by Kawabe (1978) are the same within uncertainty. Other quartz polymorphs are not stable in the T range covered by this database. In addition, no difference in oxygen isotope fractionation between various SiO 2 polymorphs has been measured in natural samples. Thus, only quartz was included in the database and its calibration can be applied to other SiO 2 polymorphs.

Feldspar group
Oxygen isotope fractionation involving albite, anorthite and K-feldspar has been widely investigated with different methods, resulting in a well-documented compilation of fractionation factors (Table 1) 4. Secondary data and internally consistent fractionation functions for couples involving selected anhydrous silicates (see text for discussion). Secondary data are reported as blue circles with vertical error bars (assigned fixed uncertainty of 0Á3&, see text), internally consistent fractionation functions as black lines and uncertainty envelopes (at 2r) as grey fields. Histograms associated with each plot show the distribution of the secondary data with respect to our model as vertical distance between each point and the model. Richter & Hoernes, 1988;Zheng, 1993a) are in agreement with each other within 1&. Data for quartz/feldspar exhibit a scatter caused by small differences in the quartz/water fractionation functions used in each study Matsuhisa et al., 1979;Matthews et al., 1983b;Chiba et al., 1989;Zheng, 1993a) (Supplementary Data 3, page 41). Variations among primary data for the couple quartz/anorthite are up to 2Á5& at T ¼ 200 C and < 1& above 600 C (Fig. 4b). The range of primary data for quartz/albite and quartz/Kfeldspar is within 0Á5&.
Our fractionation models lie within the range of the secondary data for quartz/anorthite as shown by the histogram in Fig. 4b; the primary data from Javoy et al. (1970) are outside the range of uncertainty at T < 300 C. Our model for quartz/albite lies in between a bimodal distribution of primary data (Fig. 4c), thus slightly below (0Á1-0Á3&) the calibrations of Zheng (1993a) and Chiba et al. (1989), but above (0Á3-0Á6&) the experimental data of Matsuhisa et al. (1979) and Matthews et al. (1983b) (Fig. 4c). The fractionation calculated for quartz/ K-feldspar is within the range of the primary data at T > 400 C and up to 0Á4& lower at lower T.

Garnet group
Beside semi-empirical data, natural and few experimental calibrations for garnet end-members are also available Chacko et al., 2001;Valley et al., 2003), but both approaches are limited for the purpose of this study as only data for garnet close to the composition of a pure end-member were used for fitting. The considered semi-empirical and experimental data agree with each other within 1&, with the exception of spessartine/water, for which the data disagree by 1Á5& (Fig. 4d).

Pyroxene group
Primary experimental and semi-empirical data for clinopyroxene are in good agreement (Fig. 4e), with a maximum mismatch of $ 1& between primary data for quartz/jadeite (see Supplementary Data 3, page 156). Our models for clinopyroxene are centred with respect to the primary data for the clinopyroxene/water and quartz/clinopyroxene couples; secondary data from Matthews et al. (1983a) systematically lie below the model and are partially outside the uncertainty range (see Fig. 4e).
Few data are available for fractionation involving the orthopyroxene end-members enstatite and ferrosilite (Richter & Hoernes, 1988;Zheng, 1993a) (Table 1). The divergence between primary data at low T is not significant since orthopyroxene is not stable at T < 500 C. At T > 500 C, our models are centred with respect to the data and the secondary data points are within 0Á5&.

Olivine group
Semi-empirical calculations by Richter & Hoernes (1988) and Zheng (1993a) describing the fractionation between fayalite/water and forsterite/water agree within 0Á5& between 400 and 800 C despite having an opposite concavity (Fig. 4f). Those describing the fractionation between quartz/forsterite and quartz/fayalite agree within 1& at T > 500 C. Experimental data for fractionation between quartz/forsterite and calcite/forsterite Zhang et al., 1994) are in line with semiempirical calculations for this temperature range. Our models diverge from Richter & Hoernes (1988) and Zheng (1993a) toward fractionations approaching 0& with increasing T. Our fractionation models involving tephroite are reproducing the secondary data within 0Á5& at T > 500 C.
Our models for Al 2 SiO 5 /water are close to the functions of Zheng (1993c) and Richter & Hoernes (1988). Secondary data from Hoffbauer et al. (1994) and Tennie et al. (1998) for calcite/ Al 2 SiO 5 lie systematically above our models (1-1Á5&) and part of them are out the uncertainty window. Our models for quartz/ Al 2 SiO 5 approach the data of Sharp (1995) and diverge from Zheng (1993c) up to 0Á8& with increasing T (Fig. 4g, Supplementary Data 3, pages 177-185). The calculated fractionation between andalusite/kyanite and andalusite/sillimanite are < 0Á8& at T > 200 C.

Zircon
For fractionation between quartz/zircon, natural data by Valley et al. (2003) agree within 0Á6& with the model of Zheng (1993a). The experimental study of Trail et al. (2009) is lower than the natural calibration of Valley et al. (2003) by $ 0Á5& (Fig. 4h). Semi-empirical data (Richter & Hoernes, 1988;Zheng, 1993a) are available for zircon/water fractionation (Table 1). The model of Zheng (1993a) predicts lower fractionation with a shift up to 4& at 400 C and of $ 2& at 600 C (see Supplementary Data 3, page 374).
Our model for quartz/zircon predicts a lower fractionation than Zheng (1993a), close to the data of Valley et al. (2003) (Fig. 4h). The model for zircon/water lies between the calibrations of Zheng (1993a) and Richter & Hoernes (1988).  5. Secondary data and internally consistent fractionation functions for couples involving hydrous silicates (see text for discussion). Secondary data are reported as blue circles with vertical error bars (assigned fixed uncertainty of 0Á3&, see text), internally consistent fractionation functions as black lines and uncertainty envelopes (at 2r) as grey fields. Histograms associated with each plot show the distribution of the secondary data with respect to our model as vertical distance between each point and the model.

Hydroxyl-bearing silicates Amphibole group
Available fractionation data involving amphibole are mostly limited to the semi-empirical calculations proposed by Richter & Hoernes (1988), Zheng (1993b) and Hoffbauer et al. (1994). Natural data are only available for glaucophane , and the agreement with the semi-empirical data is within 1& (Fig. 5a).
Our models are largely consistent with all these data, with a maximum difference of $ 0Á5&. The fractionation predicted by our model for quartz/glaucophane is within the range defined by semi-empirical and natural data. The model is closer to the fractionation functions of Zheng (1993b) and Javoy et al. (1970), up to 1& lower than Kohn & Valley (1998c) at T < 500 C, and follows the calibration of Javoy et al. (1970) above 500 C (Fig. 5a).

Chlorite group
Oxygen isotope fractionation factors involving chlorite group minerals have been calculated by Zheng (1993b). Natural data are available for chlorite/water (Wenner & Taylor, 1971) and quartz/chlorite fractionation (Lacroix & Vennemann, 2015), without distinction between chlorite end-members. The chlorite end-members listed above and a generic chlorite have been included in the dataset. Data for chlorite/water from Wenner & Taylor (1971) predict a significantly higher fractionation than the calculation of Richter & Hoernes (1988), up to 3& at T < 300 C (Fig. 5b).
Our models for quartz/chlorite end-members differ from the secondary data by less than 0Á5&. For chlorite/ water, our model follows the calibration of Wenner & Taylor (1971) at T < 400 C, and it is up to 0Á5& smaller than the calibration of Richter & Hoernes (1988) at higher T > 200 C (Fig. 5b). Fractionation predicted by our model for quartz/chlorite agrees within 0Á2& with the natural data of Lacroix & Vennemann (2015) (Fig. 5c).

Epidote group
Semi-empirical, experimental and natural data are available for oxygen isotope fractionation involving epidote and zoisite. The experimental calibrations by Matthews et al. (1983c) for mineral/water fractionation is up to 2& higher than the calibration of Zheng (1993b) between 500 and 800 C, but the model of Richter & Hoernes (1988) is significantly lower than the other available data below 500 C (Fig. 5d). The calibration of Zheng (1993b) for calcite/zoisite agrees with the one of  within 0Á3&, while it is $ 2Á5& higher for quartz/zoisite. Our model for zoisite/water is within the primary data range (Fig. 5d), and the one for quartz/ zoisite is close to the data of   (Supplementary Data 3, page 68). The scatter of the secondary data causes large uncertainties on the fractionation parameters (Table 1).

Mica group
The different types of primary data are mostly within 1Á5& of each other for each mica type (e.g. Fig. 5e). Notable exceptions are the fractionation functions involving margarite, with differences up to 3& (Richter & Hoernes, 1988;Zheng, 1993b;Hoffbauer et al., 1994), and paragonite/water Richter & Hoernes, 1988;Zheng, 1993b), with a mismatch up to 2& at 350 C, which decreases with increasing T.
Our model predicts intermediate fractionation values for the mica/water couples. Calculated fractionation functions for calcite/mica couples are systematically lower (0Á5-1&) than the data of Hoffbauer et al. (1994). Our model for quartz/muscovite lies in the middle of the available primary data between 500 and 800 C (Fig. 5e).

Serpentine group
A limited amount of semi-empirical (Richter & Hoernes, 1988;Zheng, 1993b) and natural data (Wenner & Taylor, 1971;Frü h-Green et al., 1996), and no experimental data for oxygen isotope fractionation are available for serpentine (Table 1). Zheng (1993b) calculated fractionation factors for amesite, lizardite and serpentine (of antigorite structure) versus water, quartz and calcite. Data for amesite/water fractionation from Zheng (1993b) are consistent within 1& with those from Richter & Hoernes (1988) at T > 400 C. Natural data for serpentine/water fractionation by Wenner & Taylor (1971) and Frü h-Green et al. (1996) agree within 1& with the calculation of Zheng (1993b) (Fig. 5f). The only dataset involving chrysotile is from Richter & Hoernes (1988) for chrysotile/water. More than one set of primary data is available for the serpentine/water and amesite/water couples, and our models lies close to the calibration of Zheng (1993b) up to 400 C, and diverge toward lower values at higher T (Fig. 5f). Our model for chrysotile/water shows an opposite concavity with respect to the primary data of Richter & Hoernes (1988), but the fit is within 0Á3& in the range of the secondary data (see Supplementary Data 3, page 338).

Calcite and dolomite
Calcite is the second most documented mineral for oxygen isotope fractionation experiments and calculations after quartz and the first among carbonates, followed by dolomite (Table 1). Data for calcite/water fractionation (Fig. 6a) show a good consistency and are mostly derived from experimental data Northrop & Clayton, 1966;. Conversely, the secondary data for quartz/calcite fractionation are more scattered (Fig. 6b), especially at low T (up to 2&). This is caused by (1) the use of independent experimental data for quartz/water and calcite/water fractionation, (2) large uncertainties on T for natural samples (e.g. Sharp & Kirschner, 1994), or (3) incomplete isotopic equilibration in experiments and natural samples (e.g. Northrop & Clayton, 1966;. The fractionation between calcite and dolomite derived from experimental and natural data (Northrop & Clayton, 1966;Sheppard & Schwarcz, 1970) is up to 3& lower (at 100 C) than the one calculated by Zheng & Bö ttcher (2016) (Fig. 6c).
Our model for calcite/water is close to the calibration of  at T > 450 C, causing a systematic shift of -1& with respect to the data of Northrop & Clayton (1966) and . Between 200 and Fig. 6. Secondary data and internally consistent fractionation functions for couples involving selected non-silicate minerals (see text for discussion). Secondary data are reported as blue circles with vertical error bars (assigned fixed uncertainty of 0Á3&, see text), internally consistent fractionation functions as black lines and uncertainty envelopes (at 2r) as grey fields. Histograms associated with each plot show the distribution of the secondary data with respect to our model as vertical distance between each point and the model. 400 C, our model approaches the calibration of   (Fig. 6a). Regarding calcite/water fractionation at low T (10-40 C), Coplen (2007) suggested that previous experimental and many biological empirical equations may suffer from a kinetic effect, underestimating the fractionation value of $ 2&. Our model in this temperature range fits the calibration of Coplen (2007) (see Supplementary Data 3, page 3), pointing to its consistency with the other data involving calcite or water. Our model for quartz/calcite has to accommodate the scattering primary data and ends up close to the fractionation function of  and Clayton et al. (1989). The secondary data of Sharp & Kirschner (1994) are above our model and partially outside the uncertainty range. By contrast, the model of Chacko et al. (1996) is significantly lower (1&) and outside the uncertainty range (Fig. 6b). Our model for calcite/dolomite is within the range of the experimental data up to 350 C, and up to 0Á8& lower at higher T (Fig. 6c). This discrepancy can be related to the abundant data available for dolomite/water fractionation that contribute to define the fractionation functions.

Iron oxides
Most of the primary data for oxygen isotope fractionation between magnetite and quartz, water or calcite agree within 2& in the range of 500 to 900 C. Only the empirical calibration of Bottinga & Javoy (1973) for magnetite/water fractionation is discordant and has been excluded from the optimization. Data for quartz/ magnetite diverge at lower T, up to 8& between the calibrations of Zheng (1991) and Chiba et al. (1989). Our model for quartz/magnetite is centred with respect to the secondary data at T > 500 C and fits the calibration of Zheng (1991) at lower T (see Supplementary Data 3, page 295). Our model for calcite/magnetite is up to 0Á5& higher than the highest batch of secondary data (Zheng, 1991), and for magnetite/water the mismatch with the data is within 1&. These deviations are caused by the abundance of data for quartz/magnetite that controls the other fractionation functions involving magnetite.
The hematite/water fractionation function at T < 200 C predicted by Zheng (1991) is dramatically lower (up to 14&) than experimental data by Yapp (1990) and Bao & Koch (1999) and the studies on natural samples by Clayton & Epstein (1961) with a scattered distribution (see Supplementary Data 3, page 283). Our model plots close to the experimental secondary data at low T, but with a large associated uncertainty derived from the scattering of the data. This large uncertainty represents a serious limitation for the modelling and use of d 18 O of hematite in geochemical studies.

Titanium oxides
Data for oxygen isotope fractionation involving rutile show a poor agreement between the semi-empirical calculations of Zheng (1991) and the more consistent sets of experimental and natural data (Addy & Garlick, 1974;Matthews & Schliestedt, 1984;Agrinier, 1991;Bird et al., 1994;Chacko et al., 2001) (i.e. quartz/rutile, Fig. 6d). In particular, semi-empirical rutile/water fractionation data (Zheng, 1991) show an opposite trend with respect to the experimental data (Addy & Garlick, 1974;, which are instead consistent with each other (see Supplementary Data 3, page 307). Considering the consistency between the calculations of Zheng (1991) and Bird et al. (1994) at low T and between the data of Addy & Garlick (1974) and  in the 450-800 C range, no primary data was excluded from the optimization.

Hydroxides
Most of the data for hydroxide minerals are from   (Table 1). New experimental and/or natural data are required for improving the quality of the present data.

Brucite
Four sets of primary data are available for brucite/water fractionation (Savin & Lee, 1988;Saccocia et al., 1998;Xu & Zheng, 1999). The results of the two semi-empirical methods from Savin & Lee (1988) and  are a remarkably different in slope and absolute vales (up to $ 5& at 200 C and $ 4& at 450 C). A detailed discussion on the equation proposed by Savin & Lee (1988) is available in . The calculation of  was preferred because of better agreement with experimental data Xu & Zheng, 1999). The dataset of Savin & Lee (1988) was, therefore, excluded for the optimization. Our model for brucite/water agrees within 1& with the data from  and Xu & Zheng (1999) at T < 150 C, while it predicts higher fractionation values (up to 1&) than  up to 600 C due to the influence of the experimental data by Saccocia et al. (1998) (Fig. 6e).

Phosphates Apatite
The only experimental data available for apatite are from Fortier & Lü ttge (1995). Chacko et al. (2001) proposed an equation obtained by linear interpolation of those experimental points. Fractionation factors involving F-apatite, Cl-apatite and OH-apatite were calculated by Zheng (1996). These different types of apatite have not been implemented in the database and the model relies on experimental data. Our model for calcite/apatite is up to 0Á5& higher than the secondary data, and the one for quartz/apatite up to 0Á3& higher (Fig. 6f,  Supplementary Data 3, page 216).

Monazite
Fractionation factors for monazite are based on the datasets by Breecker & Sharp (2007) and Rubatto et al. (2014). In the latter, the authors propose a revision to the previous calculation by Zheng (1996), which was, therefore, not considered in our calculation. Our model agrees within 0Á5& with the calibration of Breecker & Sharp (2007) at T > 500 C, being $ 0Á8& higher than the data of Rubatto et al. (2014)

Temperature range
As mentioned above, experimental or natural data at T < 200 C are available for few phases and typically exhibit a larger spread with respect to data obtained at higher T. Despite those limitations, low-T data are of key importance beside semi-empirical calculations for constraining the shape of the calculated fractionation functions, which should flatten out at increasing T because the fractionation between two phases decreases. During the global optimization, the best fit among all the data is achieved with no additional constraints on the convexity or concavity of the functions other than the secondary data. When the fractionation between two phases is small at T > 200 C, i.e. the function curvature is less pronounced, the solution that minimizes the sum of the overall residuals without considering the low-T data may result in a function with opposite curvature with respect to the primary data. Therefore, considering low-T data overall improves the match between the shape of the primary data and our model. Nevertheless, we do not recommend applying the database at T < 200 C. Furthermore, in this T range oxygen isotope fractionation between phases can be larger than 10& and the relation in Equation (5) is no longer robust. The validity of the approximation in Equation (5) and the reliability of the fit, together with the cause of possibly large uncertainties, must be evaluated by examining the results reported in Supplementary Data 3.
The choice of a maximum temperature of 900 C limits the contribution of the carbonate-exchange experiments by Clayton et al. (1989) performed at temperature up to 1000 C, those of Chiba et al. (1989), Gautason et al. (1993), , Chacko et al. (2001) and Breecker & Sharp (2007) performed up to 1200 C, and the piston-cylinder experiments of Beecker & Sharp (2007) (quartz/monazite fractionation) and Trail et al. (2009) (quartz/zircon fractionation) performed up to 1000 C. The importance of these primary data is still recognized and they are all included and discretized up to 900 C. However, they represent less than 10 % of the experimental data. Extending the T range would overemphasize the contribution of the semi-empirical calibrations at high T.

Thermometry
The temperature of equilibration between two phases i and j can be retrieved by inverting Equation (4), provided that their oxygen isotopic compositions are known. Assuming the validity of the approximation in Equation (5), T can be calculated as The inversion results in two solutions for T, of which only one is meaningful for geological interpretation (i.e. real and positive number).
A MATLABV C -based Graphical User Interface (GUI) software THERMOOX (Fig. SD4-3 in Supplementary Data 4) is distributed together with the internally consistent database. The equilibrium temperature and related uncertainties can be obtained from oxygen isotope data of two minerals among all the available ones in the database. Because the calculation performed by THERMOOX and based on our database (1) considers solid solutions, (2) determines uncertainties, (3) allows the use of phase couples for which no direct experimental or natural data are available, and (4) calculates the uncertainty also on the fractionation functions for such phase couples, it offers some advantages over existing databases compiling fractionation functions for individual mineral pairs.
The equilibrium temperature is calculated from Equation (19) and reported with three uncertainties: (1) absolute uncertainty, (2) relative uncertainty, and (3) total uncertainty (Fig. 7, Fig. SD4-3 in Supplementary Data 4). The absolute uncertainty reflects the precision of the fractionation function compared to the primary data and is determined from the uncertainties on the parameters A, B and C (Table 1). The relative uncertainty reflects the quality of the isotopic analyses and is based on the analytical uncertainties on the measured d 18 O values. The combination of the absolute and relative uncertainties results in a total uncertainty on the calculated equilibrium temperature. Because the fractionation function is not linear, the calculated upper and lower uncertainties in temperature (UAU and LAU for absolute uncertainty; URU and LRU for relative uncertainty; UTU and LTU for total uncertainty) are not symmetric.
A number of factors directly influence accuracy and precision of the results. The assumption of C ¼ 0 forces our model to diverge from primary data for which C 6 ¼ 0 (largely represented by semi-empirical data). This in turn causes an increase of the uncertainty on A and B. Additionally, the global optimization method has an effect on the different mineral/mineral fractionation factors. The multi-dimensional fit produces a database consistent with the secondary data, but it may also lead to shifts in fractionation values and larger uncertainties for specific mineral couples than otherwise obtained by fitting a single set of experiments or natural data.
The precision of oxygen isotope thermometry is related to the mathematical expression of the oxygen isotope fractionation functions, which are quadratic functions with a variation in slope (i.e. the value of the derivative) at changing T depending on the coefficients A and B. Fractionation functions having a steeper slope (D' i-j ) are less sensitive to minor changes in D 18 O iÀj , generating more precise temperature results. The effects of the slope on precision are illustrated in Fig. 8. The shift in the temperature calculated that results from a variation of 0Á2& in D 18 O iÀj at 600 C was quantified for different couples with decreasing slopes assuming the relation in Equation (5) valid. As already discussed, the fractionation function (Eq. 4) has a horizontal asymptote C ¼ 0 for T!1, so that its shape is flatter at higher T, and, therefore, the same variation in D 18 O iÀj produces larger shifts at higher T (Fig. 7). This decrease in precision affects the uncertainty of the temperature estimates for high-T mineral assemblages.
If solid solutions are involved, the molar fractions of each end-member are used to recalculate the fractionation of the mineral (see Eq. 11). The application of the fractionation model to solid solutions requires the propagation of the uncertainties on the fractionation of single end-members. Aside from this, the accuracy of the results may depend on the choice of end-member fractions. Considering the same D 18 O iÀj between quartz and different garnet end-members, the calculated equilibrium temperatures vary by a maximum of 30 C between 450 and 550 C for andradite, almandine, grossular, spessartine and uvarovite. This variation is in the same order of magnitude of the absolute Fig. 7. Uncertainty on equilibrium temperature calculated by THERMOOX. The mineral couple quartz/muscovite is taken as example. The approximation reported in Equation (5)   uncertainties of the models for each garnet endmember at the same T (20 to 30 C). Exceptions are the quartz/pyrope and quartz/melanite pairs, for which the calculated temperatures are 30 to 50 C higher and lower, respectively than the average temperature obtained by using the other end-members. Thus, for the specific case of garnet, the chosen end-member proportions have a limited effect on the result. Other mineral groups show a different behaviour. Assuming equal D 18 O iÀj between quartz/muscovite (assumed as pure Kwhite mica end-member), quartz/paragonite (assumed as pure Na-white mica end-member) and quartz/phengite (assumed as K-Mg-white mica), the calculated equilibrium T in the same temperature range for quartz/ muscovite differs by $ 40 C and $ 50 C from the result for quartz/paragonite and quartz/phengite. These differences are larger than the uncertainties of each model ($ 30 C) and the correct composition of white mica needs to be used. Such effects are even larger for the clinopyroxene group. Calculated equilibrium temperatures for quartz/acmite and quartz/jadeite are consistent with each other within uncertainty, but $ 300 C lower than those calculated using quartz/diopside or quartz/hedenbergite for the same difference in isotopic composition of each of these pairs. Significant differences are noted for plagioclase and amphiboles. In such cases, the use of solid solutions instead of pure end-members alone is of key importance to produce meaningful results.

Application of oxygen isotope fractionation factors to natural samples
The oxygen isotopic composition of a mineral in sedimentary, magmatic and metamorphic systems is controlled by several factors (e.g. Savin & Epstein, 1970): (1) the bulk isotopic composition of the reactive part of the rock (excluding detrital/inherited unreactive material, i.e. mineral relics or cores of growing porphyroblasts) and the coexisting minerals; (2) the isotopic composition of the fluid; and (3) the formation temperature or the temperature at which isotopic exchange by diffusion stopped. Therefore, the interpretation of measured values in natural samples based on modelling results that rely on our (or any other) database for oxygen isotope fractionation requires detailed petrography and petrology, including textural observation, investigation of major and trace element compositions, comparison of results from thermodynamic modelling with the observed assemblages and compositions. Insitu analysis of d 18 O at the sub-mineral scale, most commonly by ion microprobe, is a fundamental tool for oxygen isotope studies of natural samples where minerals exhibit complex textures, preserve intergranular zoning and mineral relics (e.g. Baumgartner & Valley, 2001;Russell et al., 2013;Rubatto & Angiboust, 2015;Quinn et al., 2017).
As the oxygen isotope fractionation between phases changes with T as well as the stable mineral assemblage, compositional and isotopic zoning is expected in minerals that grew over a wide T range (e.g. Kohn, 1993;Vho et al., 2020). The original isotopic zoning may be later affected by oxygen diffusion. However, constraints on diffusivity of 18 O are few and limited to some minerals (e.g. Farver & Giletti, 1985;Farver, 1989;Fortier & Giletti, 1991;Vielzeuf et al., 2005). Measured isotopic profiles within single mineral grains may help us to understand to which extent diffusion changed the original 18 O/ 16 O composition (e.g. Valley & Graham, 1993;Peck et al., 2003;Vielzeuf et al., 2005;Desbois et al., 2007;Russell et al., 2013;Higashino et al., 2019). Additionally, post-growth chemical re-equilibration of major cations does not necessarily imply oxygen isotope equilibrium.
All these factors make assessing the assumption of isotopic equilibrium in complex metamorphic rocks a challenging task. The present database can be an efficient tool to explore which mineral couples are in apparent isotopic equilibrium from the measured oxygen isotope composition of different minerals or mineral zones.
The accuracy of the fractionation functions generated by our model is directly linked to the quality of both the constraints (primary data) and our procedure for data processing. Published data reporting oxygen isotope compositions of various mineral couples from localities of known P-T conditions of equilibration can be used to test the accuracy of the fractionation factors derived from our model. Equilibrium temperatures were re-evaluated for the Yakou in Rizhao eclogite (Sulu terrane, eastern China) by combining various d 18 O mineral compositions published by Zheng et al. (2002). The calculated temperatures using our database are 425 6 20 C for quartz/muscovite, 430 6 35 C for quartz/omphacite and 460 6 20 C for quartz/garnet (both for core and rim composition, Enami et al., 1993). The results for garnet/muscovite and garnet/omphacite are consistent but associated to a large uncertainty (> 100 C) due to the limited fractionation between those mineral pairs. Temperatures obtained by different mineral couples are consistent with each other and point to a re-equilibration stage during the retrograde path as already proposed for this area (Enami et al., 1993 and references therein). A similar test was run for three samples of eclogite from Trescolmen (Adula Nappe, Swiss Alps) using the isotopic compositions of Kohn & Valley (1998b). Average temperatures of 560 6 50 C and 560 6 40 C were obtained by oxygen isotope fractionation modelling of omphacite/quartz and quartz/rutile, respectively. The results are in line with the peak metamorphic conditions of 550-650 C reported for this area (Heinrich, 1986). It is important, however, to notice that garnet and zoisite in the Trescolmen eclogite show major element zoning (Zack et al., 2002) and their d 18 O values cannot necessarily be linked to a single growth stage (see above).
To conclude, despite previous considerations regarding the possibility of dealing with large uncertainties, the internally consistent database largely matches previous calculations based on selected mineral/mineral oxygen isotope fractionation or equilibrium T obtained from independent constraints.

Limitations and alternative strategies
Calculation of an extensive, comprehensive and internally-consistent dataset is not straightforward and unavoidably relies on decisions on primary data validity and processing procedures. This section addresses potential limitations and discusses possible directions to be explored in future studies.
The problem of reliability of semi-empirical data has been discussed in the literature. Chacko et al. (1996Chacko et al. ( , 2001 showed that the increment method is most robust for anhydrous silicates, but may give inaccurate results for hydrous silicates, oxides and phosphates. Chacko & Deines (2008) argued that the increment method's predictions are also poor for carbonates. To minimize this shortfall, we adopted two strategies: (1) a relative heavier weight is given to experimental data vs. natural and semi-empirical data; and (2) the T range of our calibration has been limited to an interval that is compatible with most experimental studies, avoiding the calibration curve to be overwhelmed by semi-empirical data (see section 'Temperature range'). The final result is that mismatches between semi-empirical and other primary data are mostly within 1& also for hydrous silicates and exceptions are reported in the model results section. Incremental calculations are less reliable for oxides, and experimental and natural data are available only for a few oxides. For apatite and monazite, the available natural and experimental calibrations are integrated in our model beside semi-empirical data. However, oxides and phosphates generally occur in low modal abundance in natural rocks and thus, the uncertainty associated with their oxygen fractionation has a very minor weight when modelling multi-phase systems (see below).
A novelty in our approach is the implementation of uncertainties on isotope fractionation factors and the clarification of their statistical meaning. The results show that for the most constrained phases, a precision of 0Á5-1Á0& (2r) is achieved. This level of precision may appear large for petrological applications, but is actually comparable to the final uncertainties associated with insitu measurements by ion microprobe (i.e. obtained by adding internal and external uncertainties), which are rarely better than 0Á5& (2r), especially for minerals affected by matrix-related mass fractionation such as garnet (e.g. Page et al., 2010;Martin et al., 2014). Larger uncertainties due to discrepancies between available data must be carefully evaluated when the model is applied. Discrepancies between the model and the data that occur outside the stability field of single phases are irrelevant, as pointed out above for major minerals, such as pyroxene, feldspar and garnet.
Implementation of the presented approach and testing of alternative strategies could be considered in the future for improving the precision and the accuracy of the calculated fractionation factors.
In this study, Equation (4) was chosen as it is the most commonly used expression for describing oxygen isotope fractionation and fitting most of the experimental data (see above). It has been proposed that mineralpair fractionation at T > 300-400 C is better described by a linear relation in the space 1/T 2 (e.g. Chacko et al., 2001). However, mineral/water fractionation is not linear in the 1/T 2 space, and a linear dependence in the space 1/T better represents mineral-pair fractionation at low T (e.g. Chacko et al., 2001). Equation (4) with C ¼ 0 represents the best expression for fitting our dataset, which includes mineral/water fractionation data and covers a large temperature range. Further work could investigate the possibility of using two different linear equations, one in the space 1/T 2 at high T and one in the space 1/T for low T. This would imply choosing a T threshold valid for all the considered phases, excluding mineral/water fractionation data and might require a reassessment of published experimental data not fitted following these equations.
Further constrains could be added to the calculation, for instance, based on chemical-isotopic exchanges (e.g. , 1998b. The fractionation behaviour of some phases can be determined through combined chemical and isotopic exchanges based on the known fractionation factors of other phases. This method could, for example, provide 'constraints' to end-members for which no good experimental or semi-empirical data are available. On the other hand, this approach raises the problem of treating primary data not fully consistent with chemical-isotopic exchanges and quantifying the potential bias introduced by selecting the 'well known' phases. It has to be mentioned that deviations from chemical-isotopic exchange constraints in our model are generally within the given uncertainties. A different approach to reduce the reliance on semiempirical calibrations would be to use semi-empirical data only in the temperature ranges covered by experimental or empirical data, assigning them suitable uncertainties and then regressing data using weighted regressions. This would likely require a linear fit in the space 1/T 2 and an intercept fixed at 0, instead of Equation (4) chosen for this calculation, as a large set of sparse data points may not provide enough constraints for a correct estimate of two parameters (i.e. A and B). Additional attempts could be made by also considering data for mineral pairs beside quartz/mineral and calcite/ mineral, although they are a minority. Those data mostly come from natural samples with their associated geological uncertainty.

PERSPECTIVE: COUPLING THERMODYNAMIC AND OXYGEN FORWARD MODELLING
The database for oxygen isotope fractionation between minerals that we present can be applied in a thermodynamic framework to model the evolution of d 18 O in minerals. Thermodynamic modelling can be combined with d 18 O modelling to calculate changes in d 18 O of minerals in addition to their composition and modal abundance in metamorphic systems (Kohn, 1993;Vho et al., 2020). At fixed P, T and bulk-rock composition, the mineral assemblage, mineral modes and compositions can be determined using Gibbs free energy minimization (e.g. de Capitani & Brown, 1987). Knowledge of the mineral modes, mineral composition and T evolution allow the isotope partitioning among the phases of the system to be calculated using Equation (11) at every step of the rock evolution. The conservation of the d 18 O in the system can be described as where d 18 O sys is the isotopic composition of the system (d 18 O bulk), N sys the total number of moles of oxygen in the system, p is the number of phases, M k the number of moles of phase k, N k the its number of oxygen and d 18 O k its oxygen isotope composition (Kohn, 1993). When the differential of Equation (20) is equal to 0, the system is maintained closed with respect to oxygen isotopes. Open-system behaviour can be modelled either by setting the differential of Equation (20) (Vho et al., 2020). An example is given in Table 3. A bulk d 18 O of 12Á0& was attributed to an 'average pelite' (composition from Ague, 1991). The stable mineral assemblage, composition and modal abundances of minerals and their d 18 O values have been modelled at granulite-facies conditions (stage 1, 0Á9 GPa and 800 C) and subsequently at eclogite-facies conditions (stage 2, 2Á0 GPa and 550 C). A closed system with respect to both oxygen isotope and chemical composition has been assumed. Thermodynamic modelling was performed by using Theriak-Domino (de Capitani & Brown, 1987;de Capitani & Petrakakis, 2010) and the internally consistent dataset of Holland & Powell (1998) (tc55, distributed with Theriak-Domino 04.02.2017). No mineral fractionation was considered, i.e. the whole volume of rock was assumed to fully re-equilibrate at any change in P-T conditions. The oxygen isotope composition of each mineral phase was calculated using the fractionation factors reported in Table 1. The change in the d 18 O value of each mineral between the two stages depends on changes in T and mineral assemblage. Assuming a closed system, the largest changes between stage 1 and stage 2 are predicted for rutile ($ 1Á4&). The d 18 O value of garnet is modified by only $ 0Á8& from stage 1 to stage 2. This demonstrates that change in T alone can only produce a limited shift (i.e. 1Á0&) in the oxygen isotope composition of garnet, in agreement with published estimates (e.g. Kohn, 1993;Russell et al., 2013). Further processes (e.g. mineral and fluid fractionation, input of external fluid of different isotopic composition) are needed to produce larger variations. The predicted d 18 O values of cogenetic phases for a given bulk d 18 O can be compared with data from mineral and mineral-zone analyses in natural samples. Deviations from the models may identify open-system processes, which can then be further investigated with additional isotope and petrologic techniques (Kohn, 1993). These results provide a theoretical basis for evaluating to what extend a rock has evolved either as an open system with respect to oxygen isotopes or as a closed system fully or partially re-equilibrated.

CONCLUSIONS
Advanced applications of oxygen isotope fractionation to natural samples, such as combined thermodynamic and oxygen isotope modelling, requires consideration of the entire mineral assemblages rather than individual mineral pairs. We demonstrate that it is possible to refine the oxygen isotope fractionation factors between several minerals in a consistent way. A considerable effort was made for critical assessment of all the relevant primary data and for their transformation into an internally consistent database including 153 mineral phases and a pure water fluid phase. The main steps needed to produce the presented database are: 1. The selection of primary data (fractionation factors) from published semi-empirical, experimental and natural studies of oxygen isotope fractionation between mineral/mineral and mineral/water couples; 2. The generation of secondary data points by discretizing the fractionation function for each phase couple along the temperature range [T min , T max ] and using a different density of temperature steps T inc as a weighting factor; 3. The fitting of the secondary data and the global optimization by least-square to achieve internal consistency; 4. The calculation of uncertainties by using a Monte-Carlo technique.
Comparisons with published data on oxygen isotope compositions of minerals in select natural samples indicate that the fractionation factors calculated with this approach are robust in that they successfully match natural observations. In order to facilitate calculations, the program THERMOOX was developed to compute equilibrium temperature modelling, also considering solid solutions, and to calculate the related uncertainties.
The obtained database can be applied to model the evolution of d 18 O in minerals through their metamorphic history in the recommended T range of 200-900 C, and provides a theoretical base for evaluating to what extent a rock has evolved either as an open system with respect to oxygen isotopes or as a closed system fully or partially re-equilibrated.