MAGEMin_C.jl: Saturation models

Info

• Saturation models

• E.1 Zirconium saturation

• E.2 Zirconium and sulfur saturation models

• E.3 P2O5 saturation model example

• E.4 Saturation models and bulk correction

• E.5 CO2 saturation model example

Introduction

• When the concentration of the element in the melt increases above the saturation threshold the key phase carrying the element (see Table below) is crystallized (when the bulk-rock composition allows it) which brings back the concentration of the element in the melt to saturation.

• Although trace elements are not currently partitioned among mineral phases under subsolidus conditions, accessory minerals predicted by saturation models are still incorporated into the simulation. As a result, the model keeps track of the modal abundance and evolution of these accessory phases even below the solidus. This allows for modeling of element sequestration and mineral stability, and ensures that accessory mineral fractions are accurately represented for instance during fractional crystallization.

• By default the bulk-rock composition is not directly adjusted to account for the extra elements (e.g., Fe for sulfide). However, MAGEMin_C provides the correction to apply to the bulk-rock composition to accomodate for the extra element(s) in the trace-element output structure: Out_TE.bulk_cor_wtandOut_TE.bulk_cor_mol.

Important

While saturation is computed based on available models, if the bulk-rock composition cannot satisfy the stoichiometry of the phase, the saturation level will be increased so that the element remains in the melt.

For instance, if P2O5 concentration in the melt is greater than the saturation level but there is not enough CaO to crystallize all the excess P2O5 as fluorapatite, then the saturation level of P2O5 will be increased such that only the available CaO from the melt is used to crystallize the right fluorapatite fraction.

This approach enforces mass conservation.

• Watson & Harrison, 1983 - ZrSat_model = "WH"
• Boehnke et al., 2013 - ZrSat_model = "B"
• Crisp and Berry 2022 - ZrSat_model = "CB"

Element	phase	acronym	formula	correction	TE_prediction(; arg)
Zr	zircon	zrc	ZrSiO4	SiO2 and O	ZrSat_model = "CB", "B", "WH"
S	sulfide	sulf	FeS	FeO and O	SSat_model = "Oneill21", "Liu07"
P2O5	fluorapatite	fapt	Ca5(PO4)3F	CaO (F is omitted)	P2O5Sat_model = "Klein26"
CO2	CO2 fluid	flC	CO2	CO2 (re-enters CO2 budget)	CO2Sat_model = "SY26"

E.1 Zirconium saturation

Let's first model zirconium saturation using Crisp and Berry (2022) model. Start as usual by load MAGEMin_C, declaring the database, pressure and temperature conditions, bulk-rock composition and system unit:

using MAGEMin_C
dtb     = "mp"
data    = Initialize_MAGEMin(dtb, verbose=-1, solver=0);
P,T     = 6.0, 699.0
Xoxides = ["SiO2";  "TiO2";  "Al2O3";  "FeO";   "MnO";   "MgO";   "CaO";   "Na2O";  "K2O"; "H2O"; "O"];
X       = [58.509,  1.022,   14.858, 4.371, 0.141, 4.561, 5.912, 3.296, 2.399, 10.0, 0.2];
sys_in  = "wt"

Subsequently, we need to define the element, phase, partiton coefficient and initial concentration in ug/g:

el      = ["Zr"]
ph      = ["zrc"]
KDs     = ["0.0"]   # phase crystallized from saturation models have 0.0 KDs
C0      = [400.0]   # starting concentration of elements in ppm (ug/g)

Note

• Although Zr is not partitioned into zrc, it is necessary to declare the main Zr-bearing phase and assign a "dummy" KD value of 0. This ensures that zrc will be crystallized by stoichiometry when the Zr concentration in the liquid exceeds the saturation threshold.

Then create the partition coefficient database:

KDs_dtb = create_custom_KDs_database(el, ph, KDs)

Perform the stable phase equilibrium:

out    = single_point_minimization(P, T, data, X=X, Xoxides=Xoxides, sys_in=sys_in, name_solvus=true);

And execute the partitioning of the trace-element together with the Zrsaturation model ZrSat_model = "CB":

out_TE  = TE_prediction(out, C0, KDs_dtb, dtb; 
                        ZrSat_model     = "CB")

Finalize_MAGEMin(data)

Unfolding the out_TE structure yields:

julia> out_TE.

C0            Cliq          Cmin          Csol
Sat_P2O5_liq  Sat_S_liq     Sat_Zr_liq    bulk_D
bulk_cor_mol  bulk_cor_wt   elements      fapt_wt
liq_wt_norm   ph_TE         ph_wt_norm    sulf_wt
zrc_wt

Note

fapt_wt and sulf_wt are set to NaN as no saturation model has been provided.

Warning

The zircon fraction out_TE.zrc_wt is expressed in wt fraction (wtf) and not wt% (wt% = wtf/100.0 ...).

E.2 Zirconium and sulfur saturation models

The following example shows how to modify example E.1 to add simultaneously sulflur saturation using Liu et al. (2007) model:

using MAGEMin_C
data    = Initialize_MAGEMin("mp", verbose=-1, solver=0);
P,T     = 6.0, 699.0
Xoxides = ["SiO2";  "TiO2";  "Al2O3";  "FeO";   "MnO";   "MgO";   "CaO";   "Na2O";  "K2O"; "H2O"; "O"];
X       = [58.509,  1.022,   14.858, 4.371, 0.141, 4.561, 5.912, 3.296, 2.399, 10.0, 0.2];
sys_in  = "wt"

el      = ["Zr","S"]
ph      = ["zrc","sulf"]
KDs     = ["0.0" "0.0"; "0.0" "0.0"]   # phase crystallized from saturation models have 0.0 KDs

C0      = [400.0, 100.0]   # starting concentration of elements in ppm (ug/g)
dtb     = "mp"

KDs_dtb = create_custom_KDs_database(el, ph, KDs)

out    = single_point_minimization(P, T, data, X=X, Xoxides=Xoxides, sys_in=sys_in, name_solvus=true);
out_TE  = TE_prediction(out, C0, KDs_dtb, dtb; 
                        ZrSat_model     = "CB",
                        SSat_model      = "Liu07")

Finalize_MAGEMin(data)

E.3 P2O5 saturation model example

using MAGEMin_C
dtb     = "mp"
data    = Initialize_MAGEMin(dtb, verbose=-1, solver=0);
P,T     = 6.0, 800.0
Xoxides = ["SiO2";  "TiO2";  "Al2O3";  "FeO";   "MnO";   "MgO";   "CaO";   "Na2O";  "K2O"; "H2O"; "O"];
X       = [58.509,  1.022,   14.858, 4.371, 0.141, 4.561, 5.912, 3.296, 2.399, 10.0, 0.2];
sys_in  = "wt"

el      = ["P2O5"]
ph      = ["fapt"]  #phase bearing P2O5 is fluor-apatite
KDs     = ["0.0"]   # phase crystallized from saturation models have 0.0 KDs
C0      = [400.0]   # starting concentration of elements in ppm (ug/g)

KDs_dtb = create_custom_KDs_database(el, ph, KDs)


out    = single_point_minimization(P, T, data, X=X, Xoxides=Xoxides, sys_in=sys_in, name_solvus=true);

out_TE  = TE_prediction(out, C0, KDs_dtb, dtb; 
                        P2O5Sat_model     = "Klein26",)

Finalize_MAGEMin(data)

E.4 Saturation models and bulk correction

This example shows how to correct the bulk-rock composition in a iterative manner when accessory phase crystallized when saturation is exceeded.

using MAGEMin_C
data    = Initialize_MAGEMin("mp", verbose=-1, solver=0);
P,T     = 6.0, 699.0
Xoxides = ["SiO2";  "TiO2";  "Al2O3";  "FeO";   "MnO";   "MgO";   "CaO";   "Na2O";  "K2O"; "H2O"; "O"];
X       = [58.509,  1.022,   14.858, 4.371, 0.141, 4.561, 5.912, 3.296, 2.399, 10.0, 0.2];
X_mol, Xoxides  = convertBulk4MAGEMin(X, Xoxides,"wt","mp"); sys_in   = "mol"
X_mol ./= sum(X_mol)                                                    # normalize to 1.0

el      = ["Zr","S"]
ph      = ["zrc","sulf"]
KDs     = [ "0.0" "0.0";
            "0.0" "0.0"]                                          # phase crystallized from saturation models have 0.0 KDs

C0      = [400.0, 1000.0]                                        # starting concentration of elements in ppm (ug/g)
dtb     = "mp"

KDs_dtb = create_custom_KDs_database(el, ph, KDs)

out      = Vector{out_struct}(undef,1)
out_TE   = Vector{out_TE_struct}(undef,1)

X       = copy(X_mol)
tol     = 1e-6
res     = 1.0
n0      = 0.0
ite     = 0
while res > tol && ite < 32
    out[1]     = single_point_minimization(P, T, data, X=X, Xoxides=Xoxides, sys_in=sys_in, name_solvus=true);
    out_TE[1]  = TE_prediction(out[1] , C0, KDs_dtb, dtb; 
                            ZrSat_model     = "CB",
                            SSat_model      = "Liu07");

    X       =  X_mol .- out_TE[1].bulk_cor_mol

    res     = abs(n0 - vec_norm(out_TE[1].bulk_cor_mol))
    n0      = vec_norm(out_TE[1].bulk_cor_mol)

    ite    += 1
    if ite == 32
        @warn "Saturation model did not converge in 32 iterations, residual is $res"
    end
end

which displays:

Iteration 0: residual = 0.00011674564300540276
Iteration 1: residual = 3.909415504073306e-7

Note

• out_TE[1].bulk_cor_mol is applied to the start bulk X_mol between all iterations

• The objective of the while loop is to iterate until the residual defined as the difference of the norm of the bulk correction between two iteration decrease below the tolerance threshold tol = 1e-6

E.5 CO2 saturation model example

This example shows how to use the CO₂ saturation model of Sun & Yao (2026, "SY26"). Unlike mineral saturation phases, excess CO₂ degasses into a CO₂-bearing fluid (flC) rather than crystallizing a solid. The saturation threshold is derived from the dissolved H₂O content of the melt by numerically inverting the SY26 H₂O solubility equation, then evaluating CO₂ solubility at P_CO₂ = P − P_H₂O (binary H₂O–CO₂ fluid assumption).

using MAGEMin_C
dtb     = "mp"
data    = Initialize_MAGEMin(dtb, verbose=-1, solver=0);
P,T     = 6.0, 800.0
Xoxides = ["SiO2";  "TiO2";  "Al2O3";  "FeO";   "MnO";   "MgO";   "CaO";   "Na2O";  "K2O"; "H2O"; "O"];
X       = [58.509,  1.022,   14.858, 4.371, 0.141, 4.561, 5.912, 3.296, 2.399, 10.0, 0.2];
sys_in  = "wt"

el      = ["CO2"]
ph      = ["flC"]   # CO2 degasses into a CO2-bearing fluid
KDs     = ["0.0"]   # phase formed from saturation models have 0.0 KDs
C0      = [5000.0]  # starting CO2 concentration in ppm (ug/g)

KDs_dtb = create_custom_KDs_database(el, ph, KDs)

out    = single_point_minimization(P, T, data, X=X, Xoxides=Xoxides, sys_in=sys_in, name_solvus=true);

out_TE  = TE_prediction(out, C0, KDs_dtb, dtb;
                        CO2Sat_model = "SY26")

Finalize_MAGEMin(data)

Note

• out_TE.Sat_CO2_liq gives the CO₂ saturation concentration in the melt [ppm].

• out_TE.fl_CO2_wt gives the weight fraction of CO₂ fluid formed when saturation is exceeded.

• The saturation model requires a melt phase with dissolved H₂O. If no melt or no H₂O is present, Sat_CO2_liq and fl_CO2_wt are returned as NaN.

• Unlike zircon, sulfide and apatite, the excess CO₂ is not converted to a stoichiometrically distinct solid; it re-enters the CO₂ oxide budget in out_TE.bulk_cor_wt.

Important

The SY26 model is calibrated over 1 bar – 6 GPa, 660 – 1924 °C, and covers compositions ranging from carbonatite to rhyolite. The model assumes a binary H₂O–CO₂ fluid; mixed fluids with other species (e.g., H₂S, SO₂) are not accounted for.

References

• Watson, E. B., & Harrison, T. M. (1983). Zircon saturation revisited: temperature and composition effects in a variety of crustal magma types. earth and planetary science letters, 64(2), 295-304.

• Harrison, T. M., & Watson, E. B. (1984). The behavior of apatite during crustal anatexis: equilibrium and kinetic considerations. Geochimica et cosmochimica acta, 48(7), 1467-1477.

• Bea, F., Fershtater, G., & Corretgé, L. G. (1992). The geochemistry of phosphorus in granite rocks and the effect of aluminium. Lithos, 29(1-2), 43-56

• Bockrath, C., Ballhaus, C., & Holzheid, A. (2004). Stabilities of laurite RuS2 and monosulfide liquid solution at magmatic temperature. Chemical Geology, 208(1-4), 265-271.

• Tollari, N., Toplis, M. J., & Barnes, S. J. (2006). Predicting phosphate saturation in silicate magmas: an experimental study of the effects of melt composition and temperature. Geochimica et Cosmochimica Acta, 70(6), 1518-1536.

• Liu, Y., Samaha, N. T., & Baker, D. R. (2007). Sulfur concentration at sulfide saturation (SCSS) in magmatic silicate melts. Geochimica et Cosmochimica Acta, 71(7), 1783-1799.

• Boehnke, P., Watson, E. B., Trail, D., Harrison, T. M., & Schmitt, A. K. (2013). Zircon saturation re-revisited. Chemical Geology, 351, 324-334.

• O'Neill, H. S. C. (2021). The thermodynamic controls on sulfide saturation in silicate melts with application to ocean floor basalts. Magma redox geochemistry, 177-213.

• Crisp, L. J., & Berry, A. J. (2022). A new model for zircon saturation in silicate melts. Contributions to Mineralogy and Petrology, 177(7), 71.

• Klein, B. Z., Müntener, O., Gillespie, J., & Marxer, F. (2026). Apatite saturation revisited: new model formulations and applications to igneous rocks. Contributions to Mineralogy and Petrology, 181(3), 18.

• Sun, C., & Yao, L. (2026). A unified H₂O–CO₂ solubility–speciation model for magmatic liquids: Constraints on magma storage architecture in continental rifts. Earth and Planetary Science Letters, 689, 120115.