Add remineralization using Fe

New reaction:
    ReactionReminO2_Fe_SO4_CH4

Also add secondary redox H2 consuming reactions:
    ReactionRedoxH2_O2
    ReactionRedoxH2_SO4

src/Remin.jl

@@ -14,6 +14,18 @@ const default_reminOrgOxO2 = PB.RateStoich(
     processname="remin",
 )
 
+"Stoichiometry for organic matter oxidation by FeIIIOx (0.5 Fe2O3)
+    8H+ + 4 FeIIIOx -> 4 Fe++ + 4 H2O [+ O2]
+"
+const default_reminOrgOxFeIIIOx = PB.RateStoich(
+    PB.VarProp("reminOrgOxFeIIIOx", "mol O2eq yr-1", "oxygen consumption (-ve) by remineralization",
+        attributes=(:calc_total=>true,)),
+    ((4.0, "FeIIIOx"), (-4.0, "FeII"), (8.0, "TAlk")),
+    sms_prefix="soluteflux_", 
+    sms_suffix="",
+    processname="remin",
+)
+
 "Stoichiometry and fractionation for organic matter oxidation by SO4.
  NB: normalized to O2eq"
 const default_reminOrgOxSO4 = PB.RateStoich(
@@ -316,6 +328,108 @@ function do_remin_O2_SO4_CH4(
     return nothing
 end
 
+"""
+    ReactionReminO2_Fe_SO4_CH4
+
+Organic particulate matter remineralization (O2, FeIIIOx, SO4 oxidants, remaining Corg to CH4)
+
+FeIIIreminlimit is from Van Cappellen & Wang Am J Sci (1996) https://dx.doi.org/10.2475/ajs.296.3.197,
+100 umol Fe(OH)3 g-1 assuming dry density is 3 g/cm^3 (= 3e6 g m^-3)
+
+# Parameters
+$(PARS)
+
+# Methods and Variables for default Parameters
+$(METHODS_DO)
+"""
+Base.@kwdef mutable struct ReactionReminO2_Fe_SO4_CH4{P} <: PB.AbstractReaction
+    base::PB.ReactionBase
+
+    pars::P = PB.ParametersTuple(
+        PB.ParDouble("oxreminlimit", 8e-3, units="mol m-3", 
+            description="oxygen concentration below which use of O2 for remineralisation is inhibited"),
+
+        PB.ParDouble("FeIIIOxreminlimit", 100e-6*3e6, units="mol m-3", 
+            description="FeIII oxide concentration below which use of FeIII for remineralisation is inhibited"),
+
+        PB.ParDouble("SO4reminlimit", 1000e-3, units="mol m-3", 
+            description="sulphate concentration below which use of SO4 for remineralisation is inhibited"),
+
+        PB.ParType(PB.AbstractData, "CIsotope", PB.ScalarData,
+            external=true,
+            allowed_values=PB.IsotopeTypes,
+            description="disable / enable carbon isotopes and specify isotope type"),
+        PB.ParType(PB.AbstractData, "SIsotope", PB.ScalarData,
+            external=true,
+            allowed_values=PB.IsotopeTypes,
+            description="disable / enable sulphur isotopes and specify isotope type"),
+    )
+end
+
+function PB.register_methods!(rj::ReactionReminO2_Fe_SO4_CH4)
+
+    CIsotopeType    = rj.pars.CIsotope[]
+    SIsotopeType    = rj.pars.SIsotope[]
+    @info "register_methods! $(PB.fullname(rj)) CIsotopeType=$(CIsotopeType), SIsotopeType=$(SIsotopeType)"
+
+    (input_particulate_fluxes, output_solute_fluxes) = vars_ReminParticulateFluxes(rj, CIsotopeType, add_input_Corg_delta=true)
+
+    vars = [
+        PB.VarDep("O2_conc", "mol m-3", "O2 concentration"),
+        PB.VarDep("FeIIIOx_conc", "mol m-3", "FeIII oxide concentration"),
+        # ScalarData as we only want total, not isotopes (if any)
+        PB.VarDep("SO4_conc","mol m-3", "SO4 concentration"; attributes=(:field_data=>PB.ScalarData,)),
+        default_reminOrgOxO2.ratevartemplate,
+        default_reminOrgOxFeIIIOx.ratevartemplate,
+        default_reminOrgOxSO4.ratevartemplate,
+        default_reminOrgOxCH4.ratevartemplate,
+    ]
+
+    PB.add_method_do!(
+        rj, 
+        do_remin_O2_Fe_SO4_CH4, 
+        (input_particulate_fluxes, output_solute_fluxes, PB.VarList_namedtuple(vars)),
+    )
+
+    PB.add_method_do!(rj, default_reminOrgOxO2)
+    PB.add_method_do!(rj, default_reminOrgOxFeIIIOx)
+    PB.add_method_do!(rj, default_reminOrgOxSO4, isotope_data=SIsotopeType)
+    PB.add_method_do!(rj, default_reminOrgOxCH4, isotope_data=CIsotopeType)
+
+    PB.add_method_do_totals_default!(rj)
+    PB.add_method_initialize_zero_vars_default!(rj)
+    return nothing
+end
+
+function do_remin_O2_Fe_SO4_CH4(
+    m::PB.ReactionMethod,
+    pars,
+    (input_particulate_fluxes, output_solute_fluxes, vars),
+    cellrange::PB.AbstractCellRange,
+    deltat
+)
+
+    @inbounds for i in cellrange.indices
+        O2eq = do_ReminParticulateFluxes(i, (input_particulate_fluxes, output_solute_fluxes), Ncycle=false)
+
+        freminOrgO2     = max(vars.O2_conc[i], 0.0)/(pars.oxreminlimit[] + max(vars.O2_conc[i], 0.0))
+
+        freminOrgFeIIIOx= (1.0 - freminOrgO2)*max(vars.FeIIIOx_conc[i], 0.0)/(pars.FeIIIOxreminlimit[] + max(vars.FeIIIOx_conc[i], 0.0))
+
+        freminOrgSO4    = (1.0 - freminOrgO2 - freminOrgFeIIIOx)*max(vars.SO4_conc[i], 0.0)/(pars.SO4reminlimit[] + max(vars.SO4_conc[i], 0.0))
+
+        freminOrgCH4    = 1.0 - freminOrgO2 - freminOrgFeIIIOx - freminOrgSO4
+
+        vars.reminOrgOxO2[i] = O2eq*freminOrgO2 
+        vars.reminOrgOxSO4[i] = O2eq*freminOrgSO4  # NB: rate already normalized to O2eq yr-1
+        vars.reminOrgOxFeIIIOx[i] = O2eq*freminOrgFeIIIOx  # NB: rate already normalized to O2eq yr-1
+        # CO2 -> CH4 + 2 O2eq, so each O2eq corresponds to 0.5 CH4 released
+        # this maintains redox balance, but as no N cycle imagines nitrate and CH4 are produced ! (the least bad compromise)
+        vars.reminOrgOxCH4[i] = O2eq*freminOrgCH4  
+    end
+
+    return nothing
+end
 
 
 """

src/SecondaryRedox.jl

@@ -224,6 +224,136 @@ function PB.register_methods!(rj::ReactionRedoxCH4_SO4)
 end
 
 
+"""
+    ReactionRedoxH2_O2
+
+H2 oxidation by oxygen, producing water
+
+    H2 + 0.5 O2 -> H2O
+
+# Parameters
+$(PARS)
+
+# Methods and Variables for default Parameters
+$(METHODS_DO)
+"""
+Base.@kwdef mutable struct ReactionRedoxH2_O2{P} <: PB.AbstractReaction
+    base::PB.ReactionBase
+
+    pars::P = PB.ParametersTuple(
+        #                                          (mol l-1)-1 yr-1 / l/m^3 -> (mol m-3)-1 yr-1
+        PB.ParDouble("R_H2_O2", 1e6 , units="(mol m-3)-1 yr-1",
+            description="rate constant"),
+    )
+
+    stoich_redox_H2_O2 = PB.RateStoich(
+        PB.VarProp("redox_H2_O2", "mol H2 yr-1", "H2 consumption (+ve) by H2 oxidation by O2",
+                        attributes=(:calc_total=>true,)),
+        ((-0.5, "O2"), (-1.0, "H2")), 
+        sms_prefix="", 
+        sms_suffix="_sms",
+        processname="redox",
+    )
+
+end
+
+function PB.register_methods!(rj::ReactionRedoxH2_O2)
+
+    @info "register_methods! $(PB.fullname(rj))"
+
+    ratevar = rj.stoich_redox_H2_O2.ratevartemplate
+    concvars = [
+        PB.VarDep("O2_conc", "mol m-3", "O2 concentration"),
+        PB.VarDep("H2_conc", "mol m-3", "H2 concentration"),
+    ]
+    volume = PB.VarDep("volume", "m3", "box fluid volume")
+
+    PB.add_method_do!(
+        rj,
+        do_redox_rate, 
+        (
+            PB.VarList_single(ratevar),
+            PB.VarList_tuple(concvars),
+            PB.VarList_single(volume), 
+        );
+        p=Val(:R_H2_O2), 
+        name="do_redox_H2_O2",
+    )
+
+    PB.add_method_do!(rj, rj.stoich_redox_H2_O2)
+
+    PB.add_method_do_totals_default!(rj)
+    PB.add_method_initialize_zero_vars_default!(rj)    
+    return nothing
+end
+
+
+"""
+    ReactionRedoxH2_SO4
+
+H2 oxidation by SO4, producing sulphide
+
+     H2 + 1/4 SO4-- + 2/4 H+ -> 1/4 H2S + H2O 
+
+# Parameters
+$(PARS)
+
+# Methods and Variables for default Parameters
+$(METHODS_DO)
+"""
+Base.@kwdef mutable struct ReactionRedoxH2_SO4{P} <: PB.AbstractReaction
+    base::PB.ReactionBase
+
+    pars::P = PB.ParametersTuple(
+        #                                          (mol l-1)-1 yr-1 / l/m^3 -> (mol m-3)-1 yr-1
+        PB.ParDouble("R_H2_SO4", 1e6 , units="(mol m-3)-1 yr-1",
+            description="rate constant"),
+    )
+
+    stoich_redox_H2_SO4 = PB.RateStoich(
+        PB.VarProp("redox_H2_SO4", "mol H2 yr-1", "H2 consumption (+ve) by H2 oxidation by SO4",
+                        attributes=(:calc_total=>true,)),
+        ((-0.25, "SO4"), (-1.0, "H2"), (0.25, "H2S"), (+0.5, "TAlk")), 
+        sms_prefix="", 
+        sms_suffix="_sms",
+        processname="redox",
+    )
+
+end
+
+function PB.register_methods!(rj::ReactionRedoxH2_SO4)
+
+    @info "register_methods! $(PB.fullname(rj))"
+
+    ratevar = rj.stoich_redox_H2_SO4.ratevartemplate
+    concvars = [
+        PB.VarDep("SO4_conc", "mol m-3", "SO4 concentration"),
+        PB.VarDep("H2_conc", "mol m-3", "H2 concentration"),
+    ]
+    volume = PB.VarDep("volume", "m3", "box fluid volume")
+
+    PB.add_method_do!(
+        rj,
+        do_redox_rate, 
+        (
+            PB.VarList_single(ratevar),
+            PB.VarList_tuple(concvars),
+            PB.VarList_single(volume), 
+        );
+        p=Val(:R_H2_SO4), 
+        name="do_redox_H2_SO4",
+    )
+
+    PB.add_method_do!(rj, rj.stoich_redox_H2_SO4)
+
+    PB.add_method_do_totals_default!(rj)
+    PB.add_method_initialize_zero_vars_default!(rj)    
+    return nothing
+end
+
+
+
+# calculate generic first-order rate for two concentration defined by concvars[1:2]
 function do_redox_rate(m::PB.ReactionMethod, pars, (ratevar, concvars, volume), cellrange::PB.AbstractCellRange, deltat)
     rateparnameval = m.p
     # rateparval = getfield(pars, rateparname)[]