From bebb0305edc0730a47ab957cd8245174a18068aa Mon Sep 17 00:00:00 2001 From: Daines Date: Sun, 18 Feb 2024 15:09:16 +0000 Subject: [PATCH] Add remineralization using Fe New reaction: ReactionReminO2_Fe_SO4_CH4 Also add secondary redox H2 consuming reactions: ReactionRedoxH2_O2 ReactionRedoxH2_SO4 --- src/Remin.jl | 114 ++++++++++++++++++++++++++++++++++++ src/SecondaryRedox.jl | 130 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 244 insertions(+) diff --git a/src/Remin.jl b/src/Remin.jl index 26249c6..7284e6c 100644 --- a/src/Remin.jl +++ b/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 """ diff --git a/src/SecondaryRedox.jl b/src/SecondaryRedox.jl index 8d1cf35..8686171 100644 --- a/src/SecondaryRedox.jl +++ b/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)[]