diff --git a/README.md b/README.md index f737503..726662e 100644 --- a/README.md +++ b/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. diff --git a/input_files/env.yaml b/input_files/env.yaml index 16e4d38..545c021 100644 --- a/input_files/env.yaml +++ b/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 diff --git a/src/evo/conversions.py b/src/evo/conversions.py index 6700400..2070613 100644 --- a/src/evo/conversions.py +++ b/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 # ------------------------------------------------------------------------ diff --git a/src/evo/dgs_classes.py b/src/evo/dgs_classes.py index e79c0c6..ca0dfec 100644 --- a/src/evo/dgs_classes.py +++ b/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,7 +1085,7 @@ 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 @@ -1048,7 +1093,9 @@ def iron_fraction(self, fo2, ppa=None): 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) diff --git a/src/evo/ferric.py b/src/evo/ferric.py index 056b8be..bebe8c4 100644 --- a/src/evo/ferric.py +++ b/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"] diff --git a/src/evo/readin.py b/src/evo/readin.py index 32fd13c..81a5815 100644 --- a/src/evo/readin.py +++ b/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 diff --git a/src/evo/sat_pressure.py b/src/evo/sat_pressure.py index 1d30efe..e85e169 100644 --- a/src/evo/sat_pressure.py +++ b/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, diff --git a/src/evo/solubility_laws.py b/src/evo/solubility_laws.py index 04b88eb..840c87f 100644 --- a/src/evo/solubility_laws.py +++ b/src/evo/solubility_laws.py @@ -769,12 +769,12 @@ def oneill2020(T, melt): In: Magma Redox Chemistry. """ if not bool(melt.cm_dry): - comp = melt.Cm() + comp = melt.Cm().copy() else: - comp = melt.cm_dry + comp = melt.cm_dry.copy() if comp["fe2o3"] > 0: - comp["feo"] = comp["feo"] + comp["fe2o3"] * 0.8998 + comp["feo"] = comp["feo"] + comp["fe2o3"] * 2 comp["fe2o3"] = 0.0 comp = cnvs.mol2wt(comp) @@ -815,6 +815,67 @@ def oneill2020(T, melt): return capacity +def oneill2022(T, melt): + """ + Returns the sulfate capacity of a melt. + + Parameters + ---------- + T : float + Temperature (K) + melt : Melt class + Active instance of the Melt class + + Returns + ------- + capacity : float + Sulfate (S6+) capacity of melt as a weight fraction + + References + ---------- + O'Neill, H.S.C., Mavrogenes, J.A (2022) The sulfate capacities of silicate melts. + In: GCA. + """ + + if not bool(melt.cm_dry): + comp = melt.Cm().copy() + else: + comp = melt.cm_dry.copy() + + comp = cnvs.mol2wt(comp) + + # Converts the dry wt fraction to normalized mol fraction on single cation basis. + comp = cnvs.single_cat(comp) + + Na = comp["na2o"] + K = comp["k2o"] + Ca = comp["cao"] + Mg = comp["mgo"] + Al = comp["al2o3"] + Fe2 = comp["feo"] + Mn = comp["mno"] + + capacity = ( + np.exp( + -8.02 + + ( + 21100 + + 44000 * Na + + 18700 * Mg + + 4300 * Al + + 44200 * K + + 35600 * Ca + + 12600 * Mn + + 16500 * Fe2 + ) + / T + ) + / 1000000 + ) # Convert ppm -> wt fraction + + return capacity + + # ------------------------------------------------------------------------------ # MELT CONTENT MODEL SELECTION # ------------------------------------------------------------------------------ @@ -1024,7 +1085,7 @@ def sulfate_melt(fS2, fO2, P, T, melt, run, name="nash2019"): Active instance of the Melt class run : RunDef class Active instance of the RunDef class - name : {'nash2019'} + name : {'nash2019', 'oneill2022'} The name of the solubility law to be used Returns @@ -1038,6 +1099,20 @@ def sulfate_melt(fS2, fO2, P, T, melt, run, name="nash2019"): return ratio * (sulfide_melt(fS2, fO2, P, T, melt, name=run.SULFIDE_CAPACITY)) + elif name == "oneill2022": + C6 = oneill2022(T, melt) + + if run.SULFIDE_CAPACITY == "oneill2002": + C2 = oneill2002(fO2, P, melt) + elif run.SULFIDE_CAPACITY == "oneill2020": + C2 = oneill2020(T, melt) + + K = np.exp(55921 / T - 25.07 + 0.6465 * np.log(T)) + + ratio = (C6 / C2) * K * fO2**2 + + return ratio * (sulfide_melt(fS2, fO2, P, T, melt, name=run.SULFIDE_CAPACITY)) + def sulfide_melt(fS2, fO2, P, T, melt, name="oneill2020"): """ @@ -1288,20 +1363,26 @@ def S2_fugacity( The fugacity of S2 (bar) """ - if sulfatename == "nash2019": - # sets the Fe2/Fe3 ratio and dry melt chemistry for the sulfide capacity - melt.cm_dry = melt.iron_fraction(np.log(fO2), ppa=P * 1e5)[0] + # sets the Fe2/Fe3 ratio and dry melt chemistry for the sulfide capacity + melt.cm_dry = melt.iron_fraction(np.log(fO2), ppa=P * 1e5)[0] - # need to find the ratio of s6+/s2-, then use the amount of s2- to get mS2. + if sulfidename == "oneill2002": + C2 = oneill2002(fO2, P, melt) + + elif sulfidename == "oneill2020": + C2 = oneill2020(T, melt) + + # need to find the ratio of s6+/s2-, then use the amount of s2- to get mS2. + if sulfatename == "nash2019": ratio = nash2019(fO2, P, T, melt, run) # s6+/s2- - s2_melt = s_melt * (1 / (1 + ratio)) - if sulfidename == "oneill2002": - capacity = oneill2002(fO2, P, melt) + elif sulfatename == "oneill2022": + C6 = oneill2022(T, melt) + K = np.exp(55921 / T - 25.07 + 0.6465 * np.log(T)) + ratio = (C6 / C2) * K * fO2**2 - elif sulfidename == "oneill2020": - capacity = oneill2020(T, melt) + s2_melt = s_melt * (1 / (1 + ratio)) - fS2 = (s2_melt / capacity) ** (1 / 0.5) * fO2 + fS2 = (s2_melt / C2) ** (1 / 0.5) * fO2 - return fS2 + return fS2 diff --git a/src/evo/writeout.py b/src/evo/writeout.py index 354dee0..b8ffe83 100644 --- a/src/evo/writeout.py +++ b/src/evo/writeout.py @@ -38,6 +38,7 @@ def get_data(sys, gas, melt, P): ], "fo2": [np.exp(gas.fo2[i]) for i in range(len(P))], "F": melt.F[1:] if not sys.run.SINGLE_STEP else [melt.F[-1]], + "Fe3FeT": melt.Fe3FeT[1:] if not sys.run.SINGLE_STEP else [melt.Fe3FeT[-1]], "rho_bulk": sys.rho, "rho_melt": melt.rho_store, "Exsol_vol%": [(sys.GvF[i] * 100) for i in range(len(P))],