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 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.
The centroid of the inversion parameter ranges (black diamond, MPP) is statistically the average of the initial conditions
generated for the Nelder-Mead optimizations.
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.
Middle panels: the same results projected onto the P-misfit and
T-misfit planes.
Right panel: best central model perturbation analysis results
projected onto the P-T plane.
Because all Nelder-Mead optimizations during the perturbation analysis initiate from the best central model
coordinate {P = 9994 bar, T = 999.8 K}, the scatter of the perturbation results quantifies the imprecision of the misfit function and the inversion parameter uncertainties (\(\epsilon_{\text{misfit}}\), \(\sigma_{T,\text{perturbation}}\), and \(\sigma_{P,\text{perturbation}}\)).
These are the components of the uncertainty attributable to analytical and thermodynamic measurement error.
The misfit scatter can be viewed by rotating the perturbation analysis results into the P-misfit or T-misfit planes.
When a perturbation results in a model
with lower misfit than the best central, a large blue circle
identifies the best overall model.
If the difference between the best central and best overall
model is greater than \(\epsilon_{\text{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 filtering includes only central models whose misfit is within the imprecision of the best central model (\(\epsilon_{\text{misfit}}\)), i.e., in the range [0, 2.12e-4].
The filtered central model covariance ellipse quantifies both the uncertainty of the inversion parameters (\(\sigma_{T,\text{central}}\), \(\sigma_{P,\text{central}}\)) attributable to ill-posedness of the inverse problem and the correlations between the inversion parameters.
Right panel: Filtering has no effect on the perturbation analysis results, which are identical to those shown in the previous figure.
From Eq 5, the solution to the inverse
problem is T = 999.8 ± 21.2 K, P = 9994 ± 933 bar.
Filled symbols represent models that fit all the analytical data within its uncertainty,
the fit criterion is dependent on the confidence interval (uncertainty_multiplier).
Uncertainties and correlations for additional inversion parameters, in the present case \(\sigma_{C\text{[O2]}}\), can be evaluated by selecting the parameters from the
MC_fit_plot interface.
Refer to the caption of the previous figure for additional details.¶
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_X program
meemum(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 clinopyroxene 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 component: O2
Uncertainties:
Analytical: simulated with a reciprocal function similar to Eq 8 (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_fitoptions (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 28):
\(\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: Results of the best central model perturbation analysis.
The scatter of these results quantifies the imprecision of the misfit function (\(\epsilon_{\text{misfit}}\)) and the inversion parameter uncertainties (\(\sigma_{T,\text{perturbation}}\), \(\sigma_{P,\text{perturbation}}\)).
These are the components of the total uncertainty attributable to analytical and thermodynamic measurement error.¶
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 filtering quantifies both the uncertainty of the inversion parameters (\(\sigma_{T,\text{central}}\), \(\sigma_{P,\text{central}}\)) attributable to ill-posedness of the inverse problem and the correlations between the inversion parameters.
Right panel: Filtering has no effect on the perturbation analysis results, which are identical to those shown in the previous figure.
From Eq 5, the solution to the inverse
problem is T = 1175.5 ± 30.1 K, P = 9690 ± 579 bar.
Uncertainties and correlations for additional inversion parameters, in the present case \(\sigma_{C\text{[O2]}}\), can be evaluated by selecting the parameters from the
MC_fit_plot interface.¶
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 component: O2
Uncertainties:
Analytical ([Samuel2019], NIL16.imc)
Thermodynamic ([Holland2018], hp633ver.dat)
Modified Perple_X options (perplex_option.dat): None.
Modified
MC_fitoptions (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 28):
\(\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: Results of the best central model perturbation analysis.
The scatter of these results quantifies the imprecision of the misfit function and the inversion parameter uncertainties (\(\epsilon_{\text{misfit}}\), \(\sigma_{T,\text{perturbation}}\), and \(\sigma_{P,\text{perturbation}}\)).
These are the components of the total uncertainty attributable to analytical and thermodynamic measurement error.¶
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.
Filtering quantifies both the uncertainty of the inversion parameters (\(\sigma_{T,\text{central}}\) and \(\sigma_{P,\text{central}}\)) attributable to ill-posedness of the inverse problem and the correlations between the inversion parameters.
Right panel: Filtering has no effect on the perturbation analysis results, which are identical to those shown in the previous figure.
From Eq 5, the solution to the inverse
problem is T = 920 ± 33 K, P = 8619 ± 1164.
Uncertainties and correlations for additional inversion parameters, in the present case \(\sigma_{C\text{[O2]}}\), can be evaluated by selecting the parameters from the
MC_fit_plot interface.¶
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,
a single 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 8.
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:
component: O2
limit type: coupled
limit expression:
FeO 0 0.0625(\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{bulk} = \left[0, \frac{1}{3}\right]\), Eq 28)component: H2O
limit type: simple_mass
limit expression:
simple_mass 0 0.0166(\(\text{mass} X_{\text{H2O}} = \left[0, \frac{\text{LOI}}{1 - \text{LOI}}\right]\), LOI = 0.0163, Eq 6)
Uncertainties:
Analytical (Eq 8, VC1_bulk.imc)
Thermodynamic ([Green2016], hp62ver.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 1.63 wt %)
\(\text{bulk O}_2\) = 0.153 wt %
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{bulk}\) = 0.00647 (Eq 28)
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{Mu}\) = 0.0228 (Eq 28)
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{Bio}\) = 0.00546 (Eq 28)
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{St}\) = 0.00408 (Eq 28)
\(\text{molar} \left( \frac{\text{Fe}^{3+}}{\text{Fe}^{2+}} \right)_\text{Ilm}\) = 0.0319 (Eq 28)
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.
