Fe options

README.md

@@ -32,7 +32,7 @@ and run using
 
 Or EVo can be run directly from the terminal:
 ```
-python evo input_files/chem.yaml input_files/env.yaml --output-options input_files/output.yaml
+python -m evo input_files/chem.yaml input_files/env.yaml --output-options input_files/output.yaml
 ```
 
 The model should run and produce an output file `outputs/dgs_output_*.csv`, and a set of graphs in an 'outputs' folder, if a decompression run has been selected.

input_files/env.yaml

@@ -46,9 +46,9 @@ C_MODEL: burguisser2015
 CO_MODEL: armstrong2015
 # CH4 model, pick either ardia2013, or None (preferred for higher fO2 cases above ~ IW+1)
 CH4_MODEL: ardia2013
-# S models, pick from: sulfide capacity: oneill2002, oneill2020;  sulfate capacity: nash2019  ; SCSS: liu2007
+# S models, pick from: sulfide capacity: oneill2002, oneill2020;  sulfate capacity: nash2019 , oneill2022 ; SCSS: liu2007
 SULFIDE_CAPACITY: oneill2020
-SULFATE_CAPACITY: nash2019
+SULFATE_CAPACITY: oneill2022
 SCSS: liu2007
 # N model, pick from: libourel2003 (N2 only)
 N_MODEL: libourel2003
@@ -59,6 +59,9 @@ FO2_buffer_SET: False
 FO2_buffer: FMQ
 FO2_buffer_START: 0
 
+FE3_FET_SET: False
+FE3_FET_START: 0.1
+
 FO2_SET: False
 FO2_START: 8.3777516486972E-12
 

src/evo/conversions.py

@@ -119,7 +119,7 @@ def wt2mol(c, *args):
 
 def fo2_2F(cm, t, p, lnfo2, model="kc1991"):
     """
-    Generate the fe3/fe2 ratio based on the set fO2 and fO2 model.
+    Generate the Fe2O3/FeO ratio based on the set fO2 and fO2 model.
 
     Parameters
     ----------
@@ -148,7 +148,7 @@ def fo2_2F(cm, t, p, lnfo2, model="kc1991"):
 
 def c2fo2(cm, t, p, model):
     """
-    Calculates fO2 based on Fe3/Fe2 ratio from the melt major element
+    Calculates fO2 based on Fe2O3/FeO mole ratio from the melt major element
     composition and the fO2 model.
 
     Parameters
@@ -175,13 +175,13 @@ def c2fo2(cm, t, p, model):
 
 
 def single_cat(c):
-    """Convert dry weight fraction to single cation mol fraction. All iron as FeOt"""
+    """Convert dry weight fraction to single cation mol fraction."""
 
     mol = {}
     sm = 0
 
     for x in c:
-        if x in ["al2o3", "na2o", "k2o"]:
+        if x in ["al2o3", "na2o", "k2o", "fe2o3"]:
             mol[x] = c[x] / (cnst.m[x] * 0.5)
             sm += c[x] / (cnst.m[x] * 0.5)
         elif x in ["mgo", "sio2", "cao", "tio2", "mno", "feo"]:
@@ -338,6 +338,33 @@ def generate_fo2_buffer(sys, fo2, P, buffer_choice=None):
         return fo2_2nno(np.log10(fo2), sys.T, P, sys.run.FMQ_MODEL)
 
 
+def generate_fo2_fe3fet(sys, melt, P, Fe3_fet):
+    """
+    Returns fO2 according to the fe3/fet ratio
+
+    Parameters
+    ----------
+    sys : ThermoSystem class
+        Active instance of the ThermoSystem class
+    melt : Melt class
+        Active instance of the Melt class
+    P : float
+        Pressure (pascal)
+    F : float
+        The mFe2O3/mFeO ratio
+
+    Returns
+    -------
+    float
+        Absolute fO2
+    """
+    fe3_fe2 = Fe3_fet / (1 - Fe3_fet)
+    # F = Xfe2o3/Xfeo
+    F = fe3_fe2 / 2
+
+    return ferric.kc91_fe3(melt.Cm(), sys.T, P, F=F)
+
+
 # ------------------------------------------------------------------------
 # CONVERSIONS FOR THE SOLVER SCRIPT
 # ------------------------------------------------------------------------

src/evo/dgs_classes.py

@@ -88,6 +88,8 @@ class RunDef:
         Mineral buffer to cite fO2 relative to
     FO2_buffer_START : float
         System fO2 relative to `FO2_buffer`.
+    FE3_FET: float
+        Set the redox state with the Fe3/FeT ratio
     FO2_START : float
         Set the absolute fO2. Do not use in conjunction with `FO2_buffer_START`.
     FH2_START : float
@@ -157,12 +159,16 @@ def __init__(self):
         self.SCSS = "liu2007"
         self.N_MODEL = "libourel2003"
 
-        # initial gas fugacities
+        # initial redox
         self.FO2_buffer_SET = False
         self.FO2_buffer = "FMQ"
         self.FO2_buffer_START = 0.0  # delta buffer value
+        self.FE3_FET_SET = False
+        self.FE3_FET_START = 0.0  # ratio
         self.FO2_SET = False
         self.FO2_START = 0.0  # bar
+
+        # initial gas fugacities
         self.FH2_SET = False
         self.FH2_START = 0.24  # bar
         self.FH2O_SET = False
@@ -440,14 +446,15 @@ def norm_with_gas(self, melt):  # for iron equilibration.
         # returns the 9 major oxides as weight %, with feoT and volatiles held constant.
         oxides = cnvs.norm(C_s, (1 - cons))
 
-        self.Cw.update(
-            oxides
-        )  # adds the normalised oxide values to the total system mass frac dictionary
+        # add the normalised oxide values to the total system mass frac dictionary
+        self.Cw.update(oxides)
 
-        melt.cm_dry, F = melt.iron_fraction(
-            self.FO2
-        )  # this is now a duplicate of the one inside init cons, hopefully.
+        # this is now a duplicate of the one inside init cons, hopefully.
+        melt.cm_dry, F = melt.iron_fraction(self.FO2)
         melt.F.append(F)
+        melt.Fe3FeT.append(
+            melt.cm_dry["fe2o3"] * 2 / (melt.cm_dry["fe2o3"] * 2 + melt.cm_dry["feo"])
+        )
 
     def get_wto_tot(self, melt):
         """
@@ -547,8 +554,8 @@ def fe_save(self, melt, mo2, O2):
 
         fo2 = np.log(O2.Y * mo2 * self.P)
 
-        # returns anhydrous silicate melt comp mol fraction, Fe2/Fe3.
-        F = melt.iron_fraction(fo2)[1]
+        # returns anhydrous silicate melt comp mol fraction, Fe2O3/FeO.
+        cm_dry, F = melt.iron_fraction(fo2)
 
         ofe = (
             cnst.m["o"]
@@ -560,6 +567,7 @@ def fe_save(self, melt, mo2, O2):
         self.atomicM["o"] = self.atomicM["o_tot"] - ofe
         melt.ofe.append(ofe)
         melt.F.append(F)
+        melt.Fe3FeT.append(cm_dry["fe2o3"] * 2 / (cm_dry["fe2o3"] * 2 + cm_dry["feo"]))
 
     def mass_conservation_reset(self, melt, gas):
         """
@@ -606,7 +614,9 @@ def mass_conservation_reset(self, melt, gas):
             melt.s,
             melt.n,
             melt.F,
+            melt.Fe3FeT,
             melt.ofe,
+            melt.rho_store,
         ]
         gas_wts_f = ["H2O", "O2", "H2", "CO", "CO2", "CH4", "S2", "SO2", "H2S", "N2"]
 
@@ -925,6 +935,7 @@ def __init__(self, run, sys):
         self.s = []
         self.n = []  # stores the melt n content, independent of speciation.
         self.F = []  # stores the mfe2o3/mFeO ratio at each step.
+        self.Fe3FeT = []  # stores the fe3+/fe(total) ratio
         self.ofe = []
         self.graph_current = (
             0  # number of moles of graphite in the system at current timestep
@@ -965,9 +976,37 @@ def chem_init(self, eles, chem):
             sm += tmp_Cw[e]
             length += 1
 
-        # normalise to 100 and define initial system element mass
+        # if (
+        #     self.run.FE3_FET_SET is True
+        # ):  # PL: Rather than do this, maybe store the original values to allow the
+        # ratio to be calculated later? It's not quite accurate when calculating the
+        # saturation P.
+        #     # set correct fe2o3 and feo
+        #     eles.append("fe2o3")
+
+        #     # fe3_fe2 = 1 / (1 / self.run.FE3_FET_START - 1)
+        #     fe3_fe2 = self.run.FE3_FET_START / (1 - self.run.FE3_FET_START)
+        #     F = fe3_fe2 / 2
+
+        #     tmp_Cm = cnvs.wt2mol(tmp_Cw)
+        #     tmp_Cm["fe2o3"] = F * tmp_Cm["feo"] / (2 * F + 1)
+        #     tmp_Cm["feo"] = tmp_Cm["feo"] / (1 + 2 * F)
+
+        #     # X * molecular weight
+        #     x_mw = {k: v * cnst.m[k] for k, v in tmp_Cm.items()}
+
+        #     # wt% normalised
+        #     for e in eles:
+        #         tmp_Cw[e] = x_mw[e] * 100.0 / sum(x_mw.values())
+        # else:
+        #     # normalise to 100
+        #     for e in eles:
+        #         tmp_Cw[e] = tmp_Cw[e] * 100.0 / sm
         for e in eles:
             tmp_Cw[e] = tmp_Cw[e] * 100.0 / sm
+
+        # define initial system element mass
+        for e in eles:
             self.C_s[e] = tmp_Cw[e] / 100.0 * self.run.MASS  # actual oxide mass
 
         # add Fe+++ in if component has not been specified
@@ -1000,7 +1039,13 @@ def chem_init(self, eles, chem):
 
         # if Fe3 and Fe2 exist, then need to determine system fO2;
         # also updates the Fe mass in system.
-        if "feo" in eles and "fe2o3" in eles and self.sys.FO2 is None:
+        if self.sys.FO2 is None and self.run.FE3_FET_SET:
+            self.sys.FO2 = np.log(
+                cnvs.generate_fo2_fe3fet(
+                    self.sys, self, self.sys.Ppa, self.run.FE3_FET_START
+                )
+            )
+        elif "feo" in eles and "fe2o3" in eles and self.sys.FO2 is None:
             self.sys.FO2 = np.log(
                 cnvs.c2fo2(self.Cm(), self.sys.T, self.sys.Ppa, self.run.FO2_MODEL)
             )
@@ -1040,15 +1085,17 @@ def iron_fraction(self, fo2, ppa=None):
 
         mFeoT = composition["feo"]
 
-        mfe2o3 = F * mFeoT / (1 + 0.8998 * F)
+        mfe2o3 = F * mFeoT / (1 + 2 * F)
 
         mfeo = mfe2o3 / F
 
         composition["feo"] = mfeo
 
         composition["fe2o3"] = mfe2o3
 
-        composition = cnvs.norm(composition)
+        # PL: need to think through why you would/wouldn't normalise while holding
+        #  Fe constant...
+        composition = cnvs.norm(composition)  # , cons=["feo", "fe2o3"]
         return composition, F
 
     def Cw(self):
@@ -1385,6 +1432,11 @@ def melt_composition(self, gas, mols):
             if self.run.FE_SYSTEM is False:
                 # otherwise this is done in fe_save
                 self.F.append(F)
+                self.Fe3FeT.append(
+                    self.cm_dry["fe2o3"]
+                    * 2
+                    / (self.cm_dry["fe2o3"] * 2 + self.cm_dry["feo"])
+                )
 
             self.sulfide.append(
                 sl.sulfide_melt(
@@ -1417,6 +1469,11 @@ def melt_composition(self, gas, mols):
                 fo2 = O2.Y * gas.mO2[-1] * self.sys.P
                 self.cm_dry, F = self.iron_fraction(np.log(fo2))
                 self.F.append(F)
+                self.Fe3FeT.append(
+                    self.cm_dry["fe2o3"]
+                    * 2
+                    / (self.cm_dry["fe2o3"] * 2 + self.cm_dry["feo"])
+                )
 
             self.sulfide.append(0.0)
             self.sulfate.append(0.0)

src/evo/ferric.py

@@ -77,7 +77,7 @@ def kc91_fo2(C, T, P, lnfo2):
     return F
 
 
-def kc91_fe3(Cm, T, P):
+def kc91_fe3(Cm, T, P, F=None):
     """
     Calculates the oxygen fugacity (fO2) of a melt given where the ferric/
     ferrous ratio is known.
@@ -86,11 +86,13 @@ def kc91_fe3(Cm, T, P):
     ----------
     C : dictionary
         Major element composition of the silicate melt as mole fractions
-        Required species: Al2O3, FeO, Fe2O3, CaO,  Na2O, K2O
+        Required species: Al2O3, FeO, (Fe2O3), CaO,  Na2O, K2O
     T : float
         Temperature in degrees K
     P : float
         Pressure in pascals (Pa)
+    F : float, optional
+        If provided, the fe2o3/feo ratio to use instead of the composition
 
     Returns
     -------
@@ -121,13 +123,21 @@ def kc91_fe3(Cm, T, P):
 
     T0 = 1673.0  # K
 
+    if not F:
+        assert Cm["fe2o3"] != 0.0
+        F = Cm["fe2o3"] / Cm["feo"]
+        Cm_feo = Cm["feo"] + Cm["fe2o3"] * 0.8998
+    else:
+        assert ("fe2o3" not in Cm) or (Cm["fe2o3"] == 0.0)
+        Cm_feo = Cm["feo"]
+
     FO2 = np.exp(
         (
-            np.log(Cm["fe2o3"] / Cm["feo"])
+            np.log(F)
             - b / T
             - c
             - dal2o3 * Cm["al2o3"]
-            - dfeo * (Cm["feo"] + Cm["fe2o3"] * 0.8998)
+            - dfeo * (Cm_feo)
             - dcao * Cm["cao"]
             - dna2o * Cm["na2o"]
             - dk2o * Cm["k2o"]

src/evo/readin.py

@@ -149,7 +149,9 @@ def readin_chem(f, run, sys):
 
     # test for iron/ fO2 definition
     if "feo" in ele_names and "fe2o3" not in ele_names and run.FO2_SET is not True:
-        if run.GAS_SYS == "OH" and (run.FH2_SET):
+        if run.FE3_FET_SET is True:
+            pass
+        elif run.GAS_SYS == "OH" and (run.FH2_SET):
             pass
         elif run.FH2_SET and (run.FH2O_SET or run.WTH2O_SET):
             pass

src/evo/sat_pressure.py

@@ -241,6 +241,8 @@ def get_f(P, sys, melt, gamma):
             fo2 = np.exp(
                 cnvs.generate_fo2(sys, sys.run.FO2_buffer_START, sys.run.FO2_buffer, P)
             )
+        elif run.FE3_FET_SET:
+            fo2 = cnvs.generate_fo2_fe3fet(sys, melt, P * 1e5, run.FE3_FET_START)
 
         else:
             fo2 = cnvs.c2fo2(melt.Cm(), sys.T, P * 1e5, sys.run.FO2_MODEL)
@@ -353,6 +355,11 @@ def find_p(P, sys, melt):
             sys, sys.run.FO2_buffer_START, sys.run.FO2_buffer, sys.P
         )
 
+    elif run.FE3_FET_SET:
+        sys.FO2 = np.log(
+            cnvs.generate_fo2_fe3fet(sys, melt, sys.Ppa, run.FE3_FET_START)
+        )
+
     else:
         sys.FO2 = np.log(cnvs.c2fo2(melt.Cm(), sys.T, sys.Ppa, sys.run.FO2_MODEL))
 
@@ -1067,9 +1074,9 @@ def satp_writeout(sys, melt, gas, P, values, gamma, mols, graph_sat=False):
         or sys.run.GAS_SYS == "COHSN"
     ):
         fo2 = o2y * mO2 * P
-        melt.cm_dry, F = melt.iron_fraction(
-            np.log(fo2), ppa=P * 1e5
-        )  # Set the Fe2/Fe3 ratio and dry melt chemistry for the sulfide capacity
+        # Set the Fe2/Fe3 ratio and dry melt chemistry for the sulfide capacity
+        cm_dry, F = melt.iron_fraction(np.log(fo2), ppa=P * 1e5)
+        melt.cm_dry.update(cm_dry)
 
         # use mS2 to work back and find amount of S2- in melt
         s2_melt = (
@@ -1134,6 +1141,9 @@ def satp_writeout(sys, melt, gas, P, values, gamma, mols, graph_sat=False):
         fo2_buffer: cnvs.generate_fo2_buffer(sys, (o2y * mO2 * P), P),
         "fo2": fo2,
         "F": F,
+        "Fe3FeT": melt.cm_dry["fe2o3"]
+        * 2
+        / (melt.cm_dry["fe2o3"] * 2 + melt.cm_dry["feo"]),
         "rho_bulk": rho_bulk,
         "rho_melt": rho_melt,
         "Exsol_vol%": GvF * 100,