Problem Definition
Finding the most stable phase assemblage is a challenging constrained optimization problem including both equality and inequality constraints. One has to minimize the Gibbs energy (1) of the system while satisfying the Gibbs-Duhem (2) and mass equality constraints (3) while also satisfying the mixing-in-sites inequality constraints (4).
1. Total Gibbs Energy
The Gibbs energy of a multi-component multiphase system is given by the weighted summation of the chemical potentials of all end-members and pure phases:
where
The chemical potential of a phase is either a constant for a condensed (pure) phase:
or a function for a phase within a solution:
where
For the case of ideal mixing between the end-members, the activity coefficient is unity. The mixing of a species dissolved in a condensed phase, however, rarely behaves ideally and is typically a function of both temperature and composition (mixing-on-sites formulation).
2. Gibbs-Duhem Constraint
The Gibbs-Duhem constraint is defined as:
where
3. Mass Constraint
The mass equality constraint is defined as:
where
Minimization Approach
The Gibbs minimization approach employed in MAGEMin combines discretization of the equations of state in composition space with linear programming and extends the mass-constrained Gibbs-hyperplane rotation method to account for the mixing-on-sites that takes place in silicate mineral solid solutions. For an exhaustive description of the methodology, see Riel et al. (2022).
Local minimization of solution phases
During the linear-programming stage of the legacy solver, each solution phase is minimized locally for the current chemical potentials (finding the composition with the lowest driving force). Two minimizers are available (optimizer in MAGEMin_C.jl, --ss_solver on the command line):
NLopt (default): SLSQP in the compositional variables of the phase, with the site-fraction positivity and end-member sum written as constraints.
Nullspace minimizer: BFGS restricted to the space of compositions that satisfy the linear constraints exactly. Directions are taken in the nullspace of those constraints, so they never need to be enforced during the search. Coordinates are site fractions for the THERMOCALC databases and end-member proportions for the Stixrude & Lithgow-Bertelloni databases (there, Σp = 1 holds by construction). It allocates no memory after initialization, and falls back to NLopt for any phase it cannot handle.
With the nullspace minimizer, end-members containing an oxide that is absent from the bulk are removed before the minimization. This is what allows chemical systems with absent core oxides (see Reduced chemical systems).
Algorithm Demonstration
A simplified example of the Gibbs energy minimization approach used in MAGEMin is provided at:
This MATLAB application includes two pure phases, sillimanite and quartz, and activity-composition (a-x) relations for feldspar (pl4T, Holland et al., 2021) in a reduced Na2O-CaO-K2O-Al2O3-SiO2 (NCKAS) chemical system.
Solvus phase naming
A solvus is a boundary that defines the limit of solid solubility (miscibility) of phases described by the same solution model (or a-x model). It separates a single solid phase from a region where two (or more) distinct solid phases coexist due to compositional differences.
Properly naming demixed phases is important as it allows to differentiate multiple stable instances of the same solution model e.g., plagioclase and alkali-felspar (feldspar) or muscovite and paragonite (muscovite) etc. While for some solution models, naming (or classifying) the demixed phases is relatively straightforward (e.g., for feldspar) for other solution phases, such as amphibole, the classification rules are more complex and sometimes not fully accurate.
Igneous, igneous alkali dry (ig, igad, igd)
x = SS_vec.compVariables
- spinel:
spl
if x[3] - 0.5 > 0.0; mineral_name = "cm"
elseif x[4] - 0.5 > 0.0; mineral_name = "usp"
elseif x[2] - 0.5 > 0.0; mineral_name = "mgt"
else mineral_name = "spl"- feldspar:
fsp
if x[2] - 0.5 > 0.0; mineral_name = "afs"
else mineral_name = "pl";- muscovite:
mu
if x[4] - 0.5 > 0.0; mineral_name = "pat"
else mineral_name = "mu"- amphibole:
amp
if x[3] - 0.5 > 0.0; mineral_name = "gl"
elseif -x[3] -x[4] + 0.2 > 0.0; mineral_name = "act"
else
if x[6] < 0.1; mineral_name = "cumm"
elseif -1/2*x[4]+x[6]-x[7]-x[8]-x[2]+x[3]>0.5; mineral_name = "tr"
else mineral_name = "amp"- ilmenite:
ilm
if -x[1] + 0.5 > 0.0; mineral_name = "hem"
else mineral_name = "ilm"- nepheline:
nph
if x[2] - 0.5 > 0.0; mineral_name = "K-nph"
else mineral_name = "nph"- clinopyroxene:
cpx
if x[3] - 0.6 > 0.0; mineral_name = "pig"
elseif x[4] - 0.5 > 0.0; mineral_name = "Na-cpx"
else mineral_name = "cpx"metapelite, metabasite, ultramafic (mp, mb, um, mpe, mbe, ume)
x = SS_vec.compVariables
- feldspar:
fsp
if x[2] - 0.5 > 0.0; mineral_name = "afs"
else mineral_name = "pl";- muscovite:
mu
if x[4] - 0.5 > 0.0; mineral_name = "pat"
else mineral_name = "mu"- spinel:
sp
if x[2] - 0.5 > 0.0; mineral_name = "sp"
else mineral_name = "mt"- amphibole:
amp
if x[3] - 0.5 > 0.0; mineral_name = "gl"
elseif -x[3] -x[4] + 0.2 > 0.0; mineral_name = "act"
else
if x[6] < 0.1; mineral_name = "cumm"
elseif -1/2*x[4]+x[6]-x[7]-x[8]-x[2]+x[3]>0.5; mineral_name = "tr"
else mineral_name = "amp"- ilmenite:
ilmm
if x[1] - 0.5 > 0.0; mineral_name = "ilmm"
else mineral_name = "hemm"- ilmenite:
ilm
if 1.0 - x[1] > 0.5; mineral_name = "hem"
else mineral_name = "ilm"- diopside:
dio
if x[2] > 0.0 && x[2] <= 0.3; mineral_name = "dio"
elseif x[2] > 0.3 && x[2] <= 0.7; mineral_name = "omph"
else mineral_name = "jd"- diopside:
occm
if x[2] > 0.5; mineral_name = "sid"
elseif x[3] > 0.5; mineral_name = "ank"
elseif x[1] > 0.25 && x[3] < 0.01; mineral_name = "mag"
else mineral_name = "cc"- ortho-amphibole:
oamp
if x[2] < 0.3; mineral_name = "anth"
else mineral_name = "ged"List of mineral abbreviations (for solvus)
Names get_mineral_name can return when name_solvus = true, with the parent solution phase each one is produced from.
| Acronym | Full mineral name | Returned for |
|---|---|---|
K-nph | K-nepheline | nph |
Na-cpx | Na-rich clinopyroxene | cpx |
act | Actinolite | amp |
afs | Alkali feldspar | fsp |
amp | Amphibole | amp |
ank | Ankerite | occm |
anth | Anthophyllite | oamp |
cc | Calcite | occm |
cm | Chromite | spl |
cpx | Clinopyroxene | cpx |
cumm | Cummingtonite | amp |
dio | Clinopyroxene (dio/omphacite model) | dio |
ged | Gedrite | oamp |
gl | Glaucophane | amp |
hem | Hematite | ilm |
hemm | Hematite (Mn-bearing ilm model) | ilmm |
ilm | Ilmenite | ilm |
ilmm | Ilmenite-hematite (Mn-bearing) | ilmm |
jd | Jadeite | dio |
mag | Magnesite (in carbonate context) | occm |
mgt | Magnetite | spl |
mu | Muscovite | mu |
nph | Nepheline | nph |
omph | Omphacite (Di-Jd-Hd model; Berman/Pourteau) | dio |
pat | Paragonite-Al mica component | mu |
pig | Pigeonite | cpx |
pl | Plagioclase | fsp |
sid | Siderite | occm |
smt | Magnetite-dominated spinel | sp |
sp | Spinel | sp |
spl | Spinel | spl |
tr | Tremolite | amp |
usp | Ulvöspinel | spl |
References
Holland, T. J. B., Green, E. C. R., & Powell, R. (2021). A thermodynamic model for feldspars in KAlSi3O8-NaAlSi3O8-CaAl2Si2O8 for mineral equilibrium calculations.
Riel N., Kaus B.J.P., Green E.C.R., Berlie N., (2022) MAGEMin, an Efficient Gibbs Energy Minimizer: Application to Igneous Systems. Geochemistry, Geophysics, Geosystems 23, e2022GC010427 https://doi.org/10.1029/2022GC010427