Examples¶
This section is a guide to the MC_fit example input/output files
in the examples directory.
MC_fit problems can be broadly categorized as:
Bulk equilibrium problems constrained by a measured bulk composition.
Local equilibrium problems with no unique bulk composition.
Bulk equilibrium problems constrained by measured phase proportions.
The first category is by far the least costly. The cost of the remaining two categories rises rapidly with the number of observed phases. Since any subset of equilibrium phases is also at equilibrium, it is wise, when treating problems of these latter categories, to begin with a subset of the observed phases that unambiguously represents a relict equilibrium and then to progressively add additional phases. If the addition of a phase causes the problem to become intractable, then either the phase was not part of the relict equilibrium, or the thermodynamic model for the assemblage is flawed. The resolution of such issues is called petrology.
—
Local equilibrium problems with no unique bulk composition¶
—
bl_Chi: A synthetic garnet-clinopyroxene problem¶
Fig. 3 Left panel: bl_Chi unfiltered central model results projected onto the P-T plane.
Middle panels the same results projected onto the P-misfit and
T-misfit planes.
Right panel best central model (gold-filled circle) perturbation analysis results
projected onto the P-T plane.
The centroid of the inversion parameter ranges, statistically the average of the initial conditions
generated for the Nelder-Mead optimizations of the central model analysis, is indicated by a
black diamond (MPP).
Clustering of central model results toward the upper limit of the specified P-T range is
normally an indication that the parameter ranges should be widened.
The ranges were not widened here both because the true solution was known and it was desired to
make the problem numerically challenging by having the MPP far from the true solution.
Because all Nelder-Mead optimizations during the perturbation analysis initiate from the
best central model
coordinates {P = 9994 bar, T = 999.8 K}, the scatter of the perturbation results is
due to the propagation of analytical and thermodynamic uncertainties.
This scatter is quantifies the precision of the best central model misfit,
\(\epsilon_{misfit}\), and is reported in the legend.
The misfit scatter can be visualized by rotating the perturbation analysis results
into the P-misfit or T-misfit planes.
The misfit of perturbed models tends to be greater than the misfit of the best central model,
when the perturbations result in one or more models
with lower misfit, the best overall model is identified (blue-filled circular symbol).
When there is a significant difference between the best central and best overall
model misfit it may be of interest to investigate the perturbation that leads to the
improved fit.¶
Fig. 4 left panel Filtered central model analysis results for the bl_Chi problem
illustrated in the previous figure.
The central model results have been filtered to include only models
with misfit values within \(\pm\epsilon_{\text{misfit}}\) of the best central
model misfit, i.e., [0, 2.12e-4].
The filtered central model covariance ellipse quantifies both the uncertainty of the inversion parameters attributable to ill-posedness of the inverse problem and the correlations between the inversion parameters.
The covariance ellipse for the perturbed models (right panel) quantifies the
uncertainty of inversion parameters attributable to analytical and
thermodynamic uncertainty.
From Eq 5 and the uncertainties given in the legends of the left and right
panels, the solution to the inverse
problem is T = 999.8 ± 21.2 K, P = 9994. ± 933. bar.
Refer to the caption of the previous figure for additional details.
Filled symbols represent models that fit all the analytical data within its uncertainty,
the fit criterion is dependent on the confidence interval (uncertainty_multiplier).¶
This is a synthetic problem with a known solution designed to test MC_fit
(Slides 11-15).
The problem was created by:
Computing phase equilibria for a metabasaltic composition in the system Na2O-CaO-K2O-FeO-MgO-Al2O3-TiO2-SiO2-H2O-O2 at T = 1000 K and P = 10 kbar using the
Perple_XprogramMEEMUM(bl478.dat).The compositions of coexisting clinopyroxene and garnet compositions were extracted from the
MEEMUMforward model (bl478.prn).The excess oxygen component was stripped from the garnet and clinopyroxene compositions, making O2 an unmeasured component in the inverse problem. Additionally, the forward components not present in the garnet and clinpyroxene models (K2O, TiO2, H2O) were eliminated from the original input (bl_Chi.dat).
To make the test challenging, the ranges for pressure, temperature, and the unmeasured O2 component were set so that their centroid was far from the known solution (bl_Chi.imc).
Additional details:
The analysis of this example is discussed in Strategy.
I/O files: bl_Chi
Chemical system ([Green2016], bl_Chi.dat): Na2O-CaO-FeO-MgO-Al2O3-SiO2-O2
Saturated components: None
Unmeasured components: O2
Uncertainties:
Analytical: simulated with a reciprocal function similar to Eq 6 (bl_Chi.imc). In the current version of
MC_fit, this simulation is automatic if an uncertainty is assigned a negative value.Thermodynamic: [Green2016], (hp62ver.dat).
Modified Perple_X options (perplex_option.dat): None
Modified MC_fit options (bl_Chi_MC_fit_option.dat):
number_of_tries increased to
250to ensure adequate sampling of the solution space.
Unfiltered result (Fig. 3):
\(\text{min(misfit)}\) = 1.01e-4
\(\sigma_{T,\text{perturbation}}\) = 12.5 K
\(\sigma_{P,\text{perturbation}}\) = 462. bar
\(\epsilon_{\text{misfit}}\) = 1.11e-4
Filtered result (Fig. 4):
Best central model:
\(P\) = 9994. bar ± 933. bar
\(T\) = 999.8 K ± 21.2 K
Unmeasured components (Eq 26):
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{Cpx}\) = 0.0137
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{Grt}\) = 0.0081
Dependent variables, Oxygen fugacity and Delta QFM (Appendix C):
\(\mu_{\text{O2}}\) = -442338 J/mol
\(f_{\text{O2}}\) = -14.77 bar
\(\Delta_{\text{QFM}}\) = 0.586
Bulk composition (Appendix B):
Because this problem employs neither bulk compositional nor modal constraints, the best central model bulk composition is not an appropriate model for the synthetic source rock (bl478.prn). Nonetheless, the model bulk composition can be used to calculate phase equilibria (cf. Fig. 7) with
MEEMUMorVERTEXto verify theMC_fitprediction and to graphically assess the quality of the prediction with tools such asLinaForma([Mackay2024]),IntersecT([Nerone2025]), and/orWERAMI.
The examples/bl_Chi/variants_on_bl_Chi subdirectory includes input files for calculations that explore the use of alternative misfit functions (
LSQ,wChi), Bayes score, and the influence of the compositional weighting factor for the misfit function.
—
NIL16: A natural garnet-clinopyroxene problem¶
Fig. 5 left panel: NIL16 unfiltered central model results projected onto the P-T plane.
middle panels: the same results projected onto the P-misfit and
T-misfit planes.
right panel shows the results of the perturbation analysis
of the best central model (open gold circle), the scatter of the results quantifies the
uncertainty attributable to analytical and thermodynamic uncertainty of the misfit
function (\(\pm\epsilon_{\text{misfit}}\)) and the selected inversion parameters.¶
Fig. 6 left panel: NIL16 central model results filtered to include only models with misfit
within \(\pm\epsilon_{\text{misfit}}\) of the best central model misfit.
The filtered central model results quantify the uncertainty of the inversion parameters
attributable to ill-posedness of the inverse problem.
From Eq 5, and the uncertainties given in the legends of
the left and right panels, the solution to the inverse
problem is T = 1175.5 ± 30.1 K, P = 9690. ± 579. bar.¶
Fig. 7 Calculated P–T phase relations for the best central NIL16 model composition. This
composition is not representative of bulk rock composition (Appendix B).¶
This problem treats a natural garnet-clinopyroxene assemblage from a metagabbro, the sample, NIL16-1,
abbreviated NIL16 hereafter,
is one among a large number of samples from the Nilgiri Block, Southern India, reported and analyzed by
Samuel et al. ([Samuel2019]).
Samuel et al. concluded that the Nilgiri Block was metamorphosed at
P = 7–11 kbar and T = 923–1073 K at relatively oxidizing conditions, 2.5–3 log units above the
fayalite-magnetite-quartz (QFM) buffer.
NIL16 was selected arbitrarily from the samples presented by Samuel et al.
and, for purposes of demonstration, the equilibrium assemblage simplified to include only garnet and
clinopyroxene.
In light of these simplifications, disagreement
between the MC_fit results and those of Samuel et al. should not be construed as criticism of the
latter.
A complete analysis of the NIL16 sample with MC_fit would include the additional phases and/or bulk
compositional constraints provided by Samuel et al.
—
Additional details:
I/O files: _static/examples/NIL16
Chemical system ([Holland2018], NIL16.dat): Na2O-MgO-Al2O3-TiO2-CaO-Cr2O3-FeO-SiO2-O2
Saturated components: None
Unmeasured components: O2
Uncertainties:
Analytical ([Samuel2019], NIL16.imc)
Thermodynamic ([Holland2018], hp633ver.dat)
Modified Perple_X options (perplex_option.dat): None.
Modified MC_fit options (NIL16_MC_fit_option.dat):
number_of_tries increased to
250to ensure adequate sampling of the solution space.molar_composition_input set to
Fto allow mass fraction input compositions.relative_error set to
Fto make direct use measured errors.
Unfiltered result (Fig. 5):
\(\text{min(misfit)}\) = 2.13e-2
\(\sigma_{T,\text{perturbation}}\) = 22.4 K
\(\sigma_{P,\text{perturbation}}\) = 478. bar
\(\epsilon_{\text{misfit}}\) = 4.06e-3
Filtered result (Fig. 6):
Best central model (Try 154, NIL16.out):
\(P\) = 9690. bar ± 579. bar
\(T\) = 1175.5 K ± 18.5 K
Unmeasured components (Eq 26):
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{Cpx}\) = 0.0162
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{Grt}\) = 0.00823
Dependent variables, Oxygen fugacity and Delta QFM (Appendix C):
\(\mu_{\text{O2}}\) = -422705. J/mol
\(f_{\text{O2}}\) = -9.763 bar
\(\Delta_{\text{QFM}}\) = 2.939
This result is consistent with the contention of Samuel et al. ([Samuel2019]) that the Nilgiri Block metagabbros were metamorphosed at oxidized fugacities 2-3 log units above the QFM buffer. In contrast to the
MC_fitresult, deduced entirely from silicates, Samuel et al. estimated \(\Delta_{\text{QFM}}\) from Fe-Ti oxides.
Bulk composition (Appendix B):
Because this problem employs neither bulk compositional nor modal constraints, the best central model bulk composition (Try 154, NIL16.out) is not an appropriate model for the NIL16 rock. Nonetheless, the model bulk composition can be used to calculate phase equilibria (Fig. 7) with
MEEMUMorVERTEXto verify theMC_fitprediction and to graphically assess the quality of the prediction with tools such asLinaForma([Mackay2024]),IntersecT([Nerone2025]), and/orWERAMI.
—
Bulk equilibrium problems with bulk compositional constraints¶
—
VC1_bulk: Garnet-zone metapelite with partially known bulk composition¶
Fig. 8 left panel: VC1_bulk unfiltered central model results projected onto the P-T plane.
middle panels: the same results projected onto the P-misfit and
T-misfit planes.
right panel shows the results of the perturbation analysis
of the best central model (open gold circle), the scatter of the results quantifies the
uncertainty attributable to analytical and thermodynamic uncertainty of the misfit
function (\(\pm\epsilon_{\text{misfit}}\)) and the selected inversion parameters.¶
Fig. 9 left panel: VC1_bulk central model results filtered to include only models with misfit within \(\pm\epsilon_{\text{misfit}}\) of the best central model misfit.
The filtered central model results quantify the uncertainty of the inversion parameters attributable to ill-posedness of the inverse problem.
From Eq 5, and the uncertainties given in the legends of the left and right panels, the solution to the inverse
problem is T = 920 ± 33 K, P = 8619 ± 1164.¶
Fig. 10 Phase diagram section computed for the best central model bulk composition obtained by the VC1_bulk inversion (
VC1_bulk.out and
VC1_bulk.dat). The section was computed with VERTEX and plotted with PSSECT.¶
Fig. 11 Stability field(s) of the observed VC1_bulk garnet + biotite + ilmenite + muscovite + staurolite +
kyanite + quartz assemblage in the section shown in the previous figure. PSSECT does not discriminate
between multiple phases predicted by a single solution model (e.g., Mica + Mica). The stability fields
of the observed assemblage sensu stricto (Ti-rich Ilmenite, K-rich Mica, ± Fluid) are the two highest
temperature fields at P = 7000 – 10000 bar.¶
—
This problem treats a natural, well-equilibrated, garnet-biotite-staurolite-muscovite-kyanite-quartz-ilmenite,
metapelite from Valcolla Switzerland.
The sample (VC1) is used as an example in Andrea Galli’s Applied Geothermobarometry course at the ETH.
Except for H2O and O2 components, the bulk composition is known from XRF analysis and, excepting ilmenite,
mineral compositions are known from electron microprobe analysis.
To make the problem tractable, the ilmenite composition has been guessed to be on
Ti-rich side of ilmenite-hematite solvus.
H2O and O2 are treated as unmeasured components, SiO2 is treated as a saturated component, and the unknown
analytical uncertainties are simulated with Eq 6.
Previous analysis of the VC1 sample with MC_fit made without incorporating bulk compositional constraints
explored the influence of various simplifications of the chemical model on thermobarometric
results (Slides 19-22).
—
Additional details:
I/O files: _static/examples/VC1_bulk
Chemical system ([Green2016], VC1_bulk.dat): Na2O-MgO-Al2O3-SiO2-K2O-CaO-TiO2-MnO-FeO-H2O-O2
Saturated components: SiO2
Unmeasured components: O2, H2O
Uncertainties:
Analytical (Eq 6, VC1_bulk.imc)
Thermodynamic ([Green2016], hp633ver.dat)
Modified Perple_X options (perplex_option.dat): None.
Modified MC_fit options (VC1_MC_fit_option.dat):
number_of_tries increased to
250to ensure adequate sampling of the solution space.
Unfiltered result (Fig. 8):
\(\text{min(misfit)}\) = 0.132
\(\sigma_{T,\text{perturbation}}\) = 13.4 K
\(\sigma_{P,\text{perturbation}}\) = 717. bar
\(\epsilon_{\text{misfit}}\) = 0.858
Filtered result (Fig. 9):
Best central model (Try 175, VC1_bulk.out):
\(P\) = 8620. bar ± 1164. bar
\(T\) = 920.3 K ± 33.4 K
Unmeasured components (Try 175, VC1_bulk.out):
\(\text{bulk H}_2\text{O}\) = 1.15 wt % (LOI of 2.1 wt %)
\(\text{bulk O}_2\) = 0.153 wt %
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{bulk}\) = 0.00647 (Eq 26)
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{Mu}\) = 0.0228 (Eq 26)
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{Bio}\) = 0.00546 (Eq 26)
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{St}\) = 0.00408 (Eq 26)
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{Ilm}\) = 0.0319 (Eq 26)
Dependent variables, Oxygen fugacity and Delta QFM (Appendix C):
\(\mu_{\text{O2}}\) = -411421 J/mol (Try 175)
\(f_{\text{O2}}\) = -15.380 bar
\(\Delta_{\text{QFM}}\) = 2.1987
Bulk composition (Appendix B):
Because this problem employs bulk compositional constraints, the best central model bulk composition (Try 175, VC1_bulk.out) is an appropriate model for the VC1 rock. Consequently, the phase relations of the rock can be predicted as a function of physicochemical conditions (e.g., Fig. 10 and Fig. 11).
—
Bulk equilibrium problems with modal constraints¶
TBD.
