diff --git a/docs/src/PALEOaqchem Reactions.md b/docs/src/PALEOaqchem Reactions.md index 2975f48..f2ef62a 100644 --- a/docs/src/PALEOaqchem Reactions.md +++ b/docs/src/PALEOaqchem Reactions.md @@ -23,6 +23,17 @@ Remin.ReactionReminO2_SO4_CH4 SecondaryRedox.ReactionRedoxH2S_O2 SecondaryRedox.ReactionRedoxCH4_O2 SecondaryRedox.ReactionRedoxCH4_SO4 +SecondaryRedox.ReactionRedoxFeII_O2 +SecondaryRedox.ReactionRedoxFeIII_H2S +SecondaryRedox.ReactionRedoxFeS_O2 +SecondaryRedox.ReactionRedoxFeS2pyr_O2 +``` + +## Iron-sulphur system +```@docs +FeS.ReactionFeSaq +FeS.ReactionFeSm +FeS.ReactionPyrH2S ``` ## Carbonate chemistry diff --git a/docs/src/paleoaqchem_references.bib b/docs/src/paleoaqchem_references.bib index 7b1b85c..c71096e 100644 --- a/docs/src/paleoaqchem_references.bib +++ b/docs/src/paleoaqchem_references.bib @@ -7,6 +7,18 @@ @book{Boudreau1997 year = {1997} } +@article{Dale2015, +author = {Dale, Andrew W. and Nickelsen, L. and Scholz, F. and Hensen, C. and Oschlies, A. and Wallmann, K.}, +doi = {10.1002/2014GB005017}, +journal = {Global Biogeochemical Cycles}, +month = {may}, +number = {5}, +pages = {691--707}, +title = {{A revised global estimate of dissolved iron fluxes from marine sediments}}, +volume = {29}, +year = {2015} +} + @techreport{Lewis1998, author = {Lewis, E. and Wallace, D. W. R}, institution = {Carbon Dioxide Inf. Anal. Cent., Oak Ridge Natl. Lab., Oak Ridge, Tenn.}, @@ -24,44 +36,103 @@ @article{Ozaki2011 pages = {270--279}, publisher = {Elsevier B.V.}, title = {{Conditions required for oceanic anoxia/euxinia: Constraints from a one-dimensional ocean biogeochemical cycle model}}, -url = {http://linkinghub.elsevier.com/retrieve/pii/S0012821X11000823}, volume = {304}, year = {2011} } + +@article{rickard_solubility_2006, + title = {The solubility of {FeS}}, + volume = {70}, + issn = {00167037}, + doi = {10.1016/j.gca.2006.02.029}, + number = {23}, + journal = {Geochimica et Cosmochimica Acta}, + author = {Rickard, David}, + month = dec, + year = {2006}, + pages = {5779--5789}, +} + + +@article{rickard_chemistry_2007, + title = {Chemistry of iron sulfides.}, + volume = {107}, + issn = {0009-2665}, + doi = {10.1021/cr0503658}, + number = {2}, + journal = {Chemical reviews}, + author = {Rickard, David and Luther, George W}, + month = mar, + year = {2007}, + pmid = {17261073}, + note = {ISBN: 2920874284}, + pages = {514--62}, +} + + @article{Romaniello2010a, -author = {Romaniello, Stephen J. and Derry, Louis A}, -doi = {10.1029/2009GC002712}, -journal = {Geochemistry Geophysics Geosystems}, -month = {aug}, -number = {8}, -pages = {1--18}, -title = {{Validation of an intermediate-complexity model for simulating marine biogeochemistry under anoxic conditions in the modern Black Sea}}, -url = {http://www.agu.org/pubs/crossref/2010/2009GC002712.shtml}, -volume = {11}, -year = {2010} + author = {Romaniello, Stephen J. and Derry, Louis A}, + doi = {10.1029/2009GC002712}, + journal = {Geochemistry Geophysics Geosystems}, + month = {aug}, + number = {8}, + pages = {1--18}, + title = {{Validation of an intermediate-complexity model for simulating marine biogeochemistry under anoxic conditions in the modern Black Sea}}, + volume = {11}, + year = {2010} } @article{VanCappellen1996a, -author = {{Van Cappellen}, Philippe and Wang, Y.}, -doi = {10.2475/ajs.296.3.197}, -journal = {American Journal of Science}, -month = {mar}, -number = {3}, -pages = {197--243}, -title = {{Cycling of iron and manganese in surface sediments; a general theory for the coupled transport and reaction of carbon, oxygen, nitrogen, sulfur, iron, and manganese}}, -url = {http://www.ajsonline.org/cgi/doi/10.2475/ajs.296.3.197}, -volume = {296}, -year = {1996} + author = {{Van Cappellen}, Philippe and Wang, Y.}, + doi = {10.2475/ajs.296.3.197}, + journal = {American Journal of Science}, + month = {mar}, + number = {3}, + pages = {197--243}, + title = {{Cycling of iron and manganese in surface sediments; a general theory for the coupled transport and reaction of carbon, oxygen, nitrogen, sulfur, iron, and manganese}}, + volume = {296}, + year = {1996} } @techreport{VanHeuven2011, -author = {van Heuven, S. and Pierrot, D. and Rae, J.W.B. and Lewis, E. and Wallace, D.W.R.}, -institution = {Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, Tennessee.}, -title = {MATLAB Program Developed for CO2 System Calculations. ORNL/CDIAC-105b}, -doi = {10.3334/CDIAC/otg.CO2SYS_MATLAB_v1.1}, -url = {https://cdiac.ess-dive.lbl.gov/ftp/co2sys/CO2SYS_calc_MATLAB_v1.1/}, -year = {2011}, + author = {van Heuven, S. and Pierrot, D. and Rae, J.W.B. and Lewis, E. and Wallace, D.W.R.}, + institution = {Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, Oak Ridge, Tennessee.}, + title = {MATLAB Program Developed for CO2 System Calculations. ORNL/CDIAC-105b}, + doi = {10.3334/CDIAC/otg.CO2SYS_MATLAB_v1.1}, + url = {https://cdiac.ess-dive.lbl.gov/ftp/co2sys/CO2SYS_calc_MATLAB_v1.1/}, + year = {2011}, +} + + +@article{van_de_velde_2016, + title = {The {Influence} of {Bioturbation} on {Iron} and {Sulphur} {Cycling} in {Marine} {Sediments}: {A} {Model} {Analysis}}, + volume = {22}, + issn = {1380-6165, 1573-1421}, + shorttitle = {The {Influence} of {Bioturbation} on {Iron} and {Sulphur} {Cycling} in {Marine} {Sediments}}, + doi = {10.1007/s10498-016-9301-7}, + language = {en}, + number = {5-6}, + journal = {Aquatic Geochemistry}, + author = {Van De Velde, Sebastiaan and Meysman, Filip J. R.}, + month = dec, + year = {2016}, + keywords = {fe, sediments, sulphur}, + pages = {469--504}, +} + +@article{van_de_velde_iron_2021, + title = {Iron and sulfur cycling in the {cGENIE}.muffin {Earth} system model (v0.9.21)}, + volume = {14}, + issn = {1991-9603}, + doi = {10.5194/gmd-14-2713-2021}, + language = {en}, + number = {5}, + journal = {Geoscientific Model Development}, + author = {Van De Velde, Sebastiaan J. and Hülse, Dominik and Reinhard, Christopher T. and Ridgwell, Andy}, + month = may, + year = {2021}, + pages = {2713--2745}, } @book{Zeebe2001, diff --git a/src/FeS.jl b/src/FeS.jl new file mode 100644 index 0000000..67994c9 --- /dev/null +++ b/src/FeS.jl @@ -0,0 +1,353 @@ +module FeS + +import PALEOboxes as PB +using PALEOboxes.DocStrings + +""" + ReactionFeSaq + +Equilibrium chemistry of FeII, H2S, FeS system + +Represents two equilibrium reactions as algebraic constraints to +define primary species [Fe2++] and [HS-] in terms of: +- total [S-II] = [HS-] + [H2S] + [FeSaq] +- total [FeII] = [Fe++] + [FeSaq] +given a provided fixed value of [H+] + + H2S <--> H+ + HS- eqb. const. K_H2S + FeSaq + H+ <--> Fe++ + HS- eqb. const K2_FeSaq + + +Species concentration | PALEO variable name | PALEO variable type +----------------------------|-------------------------|--------------------- +[H+] | Hp_conc | Dependency +total [S-II] | SmIIaqtot_conc | Dependency +[HS-] | Hsm_conc | State +constraint on total [S-II] | SmIIaqtot_constraint | Constraint +[H2S] | H2Ssp_conc | Property +[HS-] + [H2S] | H2Stot_conc | Property +total [FeII] | FeIIaqtot_conc | Dependency +[Fe2++] | FeII_conc | State +constraint on total [FeII] | FeIIaqtot_constraint | Constraint + + +(usually SmIIaqtot and FeIIaqtot are ODE variables and should be defined by ReactionReservoir) + +See: +- Rickard (2006) GCA https://10.1016/j.gca.2006.02.029 +- Rickard and Luther (2007) Chemical Reviews https://dx.doi.org/10.1021/cr0503658 +- van de Velde etal (2021) GMD https://10.5194/gmd-14-2713-2021 + +# Parameters +$(PARS) + +# Methods and Variables for default Parameters +$(METHODS_DO) +""" +Base.@kwdef mutable struct ReactionFeSaq{P} <: PB.AbstractReaction + base::PB.ReactionBase + + pars::P = PB.ParametersTuple( + + PB.ParDouble("K_H2S", 2.15e-7 , units="mol kg-1", + description="K_H2S = [H+][HS-]/[H2S] units mol kg-1 eqb const for H2S <--> HS- + H+ (Hofmann (2010) Aquatic Geochemistry https://10.1007/s10498-009-9084-1)"), + PB.ParDouble("K2_FeSaq", 10.0^2.2 , units="", + description="K2_FeSaq = {Fe++}[HS-]/([FeSaq][H+]) eqb constant for FeSaq + H+ <--> Fe++ + HS- (Rickard (2006) eq 8)"), + PB.ParDouble("activity_Fe2p", 0.23, units="", + description="Fe++ activity (eg van de Velde 2012 Table 3 from Davies eqn)") + + ) + + +end + + +# K_H2S from Hofmann (2010) Aquatic Geochemistry https://10.1007/s10498-009-9084-1 +# Sal = 35 +# TempK = 273.15 + 15 +# lnKH2S= (225.838 + 0.3449*sqrt(Sal) - 0.0274*Sal) - 13275.3/TempK - 34.6435*log(TempK) +# -15.35415072896447 +# exp(lnKH2S) +# 2.146728213653719e-7 # mol kg-1 + +function PB.register_methods!(rj::ReactionFeSaq) + + @info "register_methods! $(PB.fullname(rj))" + + setup_vars = [ + PB.VarState("HSm_conc", "mol m-3", "aqueous HS- concentration"), + PB.VarState("FeII_conc", "mol m-3", "aqueous Fe++ concentration"), + PB.VarConstraint("SmIIaqtot_constraint", "mol m-3", "algebraic constraint on SmIIaqtot_conc (= 0)"), + PB.VarConstraint("FeIIaqtot_constraint", "mol m-3", "algebraic constraint on FeIIaqtot_conc (= 0)"), + ] + + vars_FeSaq_eqb = vcat( + setup_vars, + [ + PB.VarDep("Hp_conc", "mol m-3", "aqueous H+ concentration"), + PB.VarDep("rho_ref", "kg m-3", "reference density"), + PB.VarDep("SmIIaqtot_conc", "mol m-3", "total aqueous S-II concentration"), + PB.VarDep("FeIIaqtot_conc", "mol m-3", "total aqueous FeII concentration"), + PB.VarProp("H2Ssp_conc", "mol m-3", "aqueous H2S species concentration"), + PB.VarProp("H2Stot_conc", "mol m-3", "aqueous H2S + HS- concentration"), + PB.VarProp("FeSaq_conc", "mol m-3", "aqueous FeS concentration"), + ] + ) + + PB.add_method_setup!( + rj, + setup_FeSaq_eqb, + (PB.VarList_namedtuple(setup_vars), ), + ) + + PB.add_method_do!( + rj, + do_FeSaq_eqb, + ( + PB.VarList_namedtuple(vars_FeSaq_eqb), + ); + ) + + return nothing +end + +function setup_FeSaq_eqb( + m::PB.ReactionMethod, + pars, + (setup_vars, ), + cellrange::PB.AbstractCellRange, + attribute_name +) + attribute_name in (:initial_value, :norm_value) || return + + norm_value = 1e-3 # for solver scaling only + + # constraint variables just need norm_value setting + if attribute_name == :norm_value + setup_vars.SmIIaqtot_constraint[cellrange.indices] .= norm_value + @info "ReactionFeSaq initialize_state:$(rpad(attribute_name,20)) $(rpad("SmIIaqtot_constraint", 30)) = $norm_value" + setup_vars.FeIIaqtot_constraint[cellrange.indices] .= norm_value + @info "ReactionFeSaq initialize_state:$(rpad(attribute_name,20)) $(rpad("FeIIaqtot_constraint", 30)) = $norm_value" + end + + # state variables need both norm_value and initial_value setting + initialnorm_value = Dict(:initial_value=>1e-9, :norm_value=>norm_value)[attribute_name] + setup_vars.HSm_conc[cellrange.indices] .= initialnorm_value + @info "ReactionFeSaq initialize_state:$(rpad(attribute_name,20)) $(rpad("HSm_conc", 30)) = $initialnorm_value" + setup_vars.FeII_conc[cellrange.indices] .= initialnorm_value + @info "ReactionFeSaq initialize_state:$(rpad(attribute_name,20)) $(rpad("FeII_conc", 30)) = $initialnorm_value" + + return nothing +end + +function do_FeSaq_eqb(m::PB.ReactionMethod, pars, (vars,), cellrange::PB.AbstractCellRange, deltat) + + K_H2S = pars.K_H2S[] # K_H2S = [H+][HS-]/[H2S] units mol kg-1 + K2_FeSaq = pars.K2_FeSaq[] + activity_Fe2p = pars.activity_Fe2p[] + + for i in cellrange.indices + # H2S <--> HS- + H+ eqb + # mol m-3 = mol m-3 * mol m-3 / (mol kg-1 * kg m-3) + vars.H2Ssp_conc[i] = vars.HSm_conc[i] * vars.Hp_conc[i] / (K_H2S * vars.rho_ref[i]) + vars.H2Stot_conc[i] = vars.HSm_conc[i] + vars.H2Ssp_conc[i] + + # FeSaq + H+ <--> Fe++ + HS- eqb + vars.FeSaq_conc[i] = (activity_Fe2p * vars.FeII_conc[i] * vars.HSm_conc[i] / vars.Hp_conc[i]) / K2_FeSaq + + # H2S budget constraint for DAE solver + vars.SmIIaqtot_constraint[i] = vars.SmIIaqtot_conc[i] - (vars.H2Stot_conc[i] + vars.FeSaq_conc[i]) + + # FeII budget constraint for DAE solver + vars.FeIIaqtot_constraint[i] = vars.FeIIaqtot_conc[i] - (vars.FeII_conc[i] + vars.FeSaq_conc[i]) + end + + return nothing +end + + +""" + ReactionFeSm + +FeS precipitation/dissolution + +Represents + + [FeSaq] <--> [FeSm] + +as a fast precipitation-dissolution reaction, using the formulation from Van Cappellen & Wang (1996) + +See: +- Van Cappellen & Wang (1996) AJS https://dx.doi.org/10.2475/ajs.296.3.197 +- Rickard (2006) GCA https://dx.doi.org/10.1016/j.gca.2006.02.029 +- Rickard and Luther (2007) Chemical Reviews https://dx.doi.org/10.1021/cr0503658 + +# Parameters +$(PARS) + +# Methods and Variables for default Parameters +$(METHODS_DO) +""" +Base.@kwdef mutable struct ReactionFeSm{P} <: PB.AbstractReaction + base::PB.ReactionBase + + pars::P = PB.ParametersTuple( + PB.ParDouble("K0_FeSm", 10.0^-5.7 , units="mol l-1", + description="solubility constant for FeSm"), + PB.ParDouble("rate_FeS_prec", 40.0 , units="mol m-3 yr-1", + description="rate of FeS precipitation"), + PB.ParDouble("rate_FeS_diss", 1e3 , units="yr-1", + description="rate of FeS dissolution"), + ) + + stoich_FeS_precdiss = PB.RateStoich( + PB.VarProp("rate_FeSm", "mol yr-1", "FeSm precipitation - dissolution rate"; + attributes=(:calc_total=>true,)), + ((-1.0, "SmIIaqtot"), (-1.0, "FeIIaqtot"), (1.0, "FeSm"), ), + sms_prefix="", + sms_suffix="_sms", + processname="mineral", + ) +end + + +function PB.register_methods!(rj::ReactionFeSm) + + @info "register_methods! $(PB.fullname(rj))" + + vars = [ + rj.stoich_FeS_precdiss.ratevartemplate, + PB.VarDep("volume_solute", "m3", "solute fluid volume"), + PB.VarDep("volume_solid", "m3", "solid-phase volume"), + PB.VarDep("FeSaq_conc", "mol m-3", "aqueous FeS concentration"), + PB.VarProp("Omega_FeSaq", "", "aqueous FeS saturation"), + PB.VarDep("FeSm_conc", "mol m-3", "FeS solid phase concentration"), + ] + + PB.add_method_do!( + rj, + do_FeSm, + ( + PB.VarList_namedtuple(vars), + ); + ) + + PB.add_method_do!(rj, rj.stoich_FeS_precdiss) + + PB.add_method_do_totals_default!(rj) + PB.add_method_initialize_zero_vars_default!(rj) + + return nothing +end + + +function do_FeSm(m::PB.ReactionMethod, pars, (vars,), cellrange::PB.AbstractCellRange, deltat) + # FeSm saturation and precipitation / dissolution + + K0_FeSm = pars.K0_FeSm[] + + for i in cellrange.indices + # FeS saturation + # mol m-3 = mol l-1 * l m-3 + FeS_sat_conc = K0_FeSm * 1000 + vars.Omega_FeSaq[i] = vars.FeSaq_conc[i] / FeS_sat_conc + + # FeSm precipitation and dissolution + # NB: bodge so AD sparsity detection sees a dependency on both FeSaq_conc and FeSm_conc + # set rate to zero first with a fake dependency on both vars, then add/subtract to set value + vars.rate_FeSm[i] = 0*vars.FeSm_conc[i] + 0*vars.FeSaq_conc[i] + + if vars.Omega_FeSaq[i] >= 1.0 # FeSm precipitation + # Van Cappellen & Wang (1996) formulation + # mol yr-1 = mol m-3 yr-1 * m3 + vars.rate_FeSm[i] += pars.rate_FeS_prec[]*(vars.Omega_FeSaq[i] - 1) * vars.volume_solute[i] + else # FeSm dissolution + # mol yr-1 = yr-1 * mol m-3 * m3 + # Van Cappellen & Wang (1996) formulation + vars.rate_FeSm[i] -= pars.rate_FeS_diss[]*(1 - vars.Omega_FeSaq[i])*max(vars.FeSm_conc[i] - 1e-12, 0.0)*vars.volume_solid[i] # limit to minimum conc of 1e-12 + end + + end + + return nothing +end + + +""" + ReactionPyrH2S + +Pyrite formation by the "H2S" (Berzelius) mechanism + + FeSm + H2S -> FeS2pyr + H2 + +with rate + + R_Pyr_H2S * [FeSm] * [H2S] + +Where R_Pyr_H2S is highly uncertain: + +[Dale2015](@cite) use R_Pyr_H2S = 1e5 M-1 yr-1 (= 1e5 *1e-3 = 100 (mol m-3-1) yr-1) + +[van_de_velde_iron_2021](@cite) use R_Pyr_H2S = 0.3708 M-1 h-1 (= 0.3708 * 8.766 (mol m-3)-1 yr-1), +where [H2S] is taken to be total H2S (HS- + H2S) + +# Parameters +$(PARS) + +# Methods and Variables for default Parameters +$(METHODS_DO) +""" +Base.@kwdef mutable struct ReactionPyrH2S{P} <: PB.AbstractReaction + base::PB.ReactionBase + + pars::P = PB.ParametersTuple( + PB.ParDouble("R_Pyr_H2S", 0.3708e0 * 8.766 , units="(mol m-3)-1 yr-1", + description="rate constant"), + ) + + stoich_FeS2pyr = PB.RateStoich( + PB.VarProp("rate_FeS2pyr_H2S", "mol yr-1", "pyrite formation rate via the H2S mechanism"; + attributes=(:calc_total=>true,)), + ((-1.0, "FeSm"), (-1.0, "SmIIaqtot"), (1.0, "FeS2pyr"), (1.0, "H2"), ), + sms_prefix="", + sms_suffix="_sms", + processname="mineral", + ) +end + +function PB.register_methods!(rj::ReactionPyrH2S) + + @info "register_methods! $(PB.fullname(rj))" + + vars = [ + rj.stoich_FeS2pyr.ratevartemplate, + PB.VarDep("FeSm_conc", "mol m-3", "FeSm concentration"), + PB.VarDep("H2Ssp_conc", "mol m-3", "H2S (species) concentration"), + PB.VarDep("volume", "m3", "box fluid volume"), + ] + + PB.add_method_do!(rj, do_PyrH2S, (PB.VarList_namedtuple(vars), )) + + PB.add_method_do!(rj, rj.stoich_FeS2pyr) + + PB.add_method_do_totals_default!(rj) + PB.add_method_initialize_zero_vars_default!(rj) + return nothing +end + + +function do_PyrH2S(m::PB.ReactionMethod, pars, (vars, ), cellrange::PB.AbstractCellRange, deltat) + + for i in cellrange.indices + # mol yr-1 + vars.rate_FeS2pyr_H2S[i] = ( pars.R_Pyr_H2S[] # (mol m-3)-1 yr-1 + * max(PB.get_total(vars.FeSm_conc[i]), 0.0) # mol m-3 + * max(PB.get_total(vars.H2Ssp_conc[i]), 0.0) # mol m-3 + * vars.volume[i]) # m3 + end + + return nothing +end + + +end diff --git a/src/PALEOaqchem.jl b/src/PALEOaqchem.jl index 3bbf38e..db334e4 100644 --- a/src/PALEOaqchem.jl +++ b/src/PALEOaqchem.jl @@ -42,6 +42,8 @@ include("Remin.jl") include("SecondaryRedox.jl") +include("FeS.jl") + include("Particle.jl") include("Boron.jl") diff --git a/src/SecondaryRedox.jl b/src/SecondaryRedox.jl index 8686171..1a31d77 100644 --- a/src/SecondaryRedox.jl +++ b/src/SecondaryRedox.jl @@ -352,6 +352,283 @@ function PB.register_methods!(rj::ReactionRedoxH2_SO4) end +################################################################################################################## +# Iron +############################################################################################################### + +""" + ReactionRedoxFeII_O2 + +FeII oxidation by oxygen, producing generic FeIIIOx (0.5 Fe2O3) + + 4 Fe++ + O2 + 4 H2O -> 4 FeIIIOx + 8 H+ + +# Parameters +$(PARS) + +# Methods and Variables for default Parameters +$(METHODS_DO) +""" +Base.@kwdef mutable struct ReactionRedoxFeII_O2{P} <: PB.AbstractReaction + base::PB.ReactionBase + + pars::P = PB.ParametersTuple( + # + PB.ParDouble("R_FeII_O2", 3.65e6/1e3 , units="(mol m-3)-1 yr-1", + description="rate constant"), + ) + + stoich_redox_FeII_O2 = PB.RateStoich( + PB.VarProp("redox_FeII_O2", "mol O2 yr-1", "oxygen consumption (+ve) by Fe++ oxidation", + attributes=(:calc_total=>true,)), + ((-1.0, "O2"), (-4.0, "FeII"), (4.0, "FeIIIOx"), (-8.0, "TAlk")), + sms_prefix="", + sms_suffix="_sms", + processname="redox", + ) + +end + +function PB.register_methods!(rj::ReactionRedoxFeII_O2) + + @info "register_methods! $(PB.fullname(rj))" + + ratevar = rj.stoich_redox_FeII_O2.ratevartemplate + concvars = [ + PB.VarDep("O2_conc", "mol m-3", "O2 concentration"), + PB.VarDep("FeII_conc", "mol m-3", "FeII 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_FeII_O2), + name="do_redox_FeII_O2", + ) + + PB.add_method_do!(rj, rj.stoich_redox_FeII_O2) + + PB.add_method_do_totals_default!(rj) + PB.add_method_initialize_zero_vars_default!(rj) + return nothing +end + + +""" + ReactionRedoxFeIII_H2S + +Sulphate-mediated iron reduction (of generic iron oxide FeIIIOx ~ 0.5 Fe2O3) + + H2S + 8 FeIIIOx + 14 H+ -> 8 Fe++ + SO4-- + 8 H2O + +NB: rate constant k [units (mol m-3)^-0.5 yr-1] = 277.2 * k [units (mol l-1)^-0.5 hr-1] + +where van de Velde (2021) give k_SMIs = 1.98e0 (mol l-1)^-0.5 hr-1 for solid-phase FeIII + +# Parameters +$(PARS) + +# Methods and Variables for default Parameters +$(METHODS_DO) +""" +Base.@kwdef mutable struct ReactionRedoxFeIII_H2S{P} <: PB.AbstractReaction + base::PB.ReactionBase + + pars::P = PB.ParametersTuple( + + PB.ParDouble("R_FeIII_H2S", 277.2*1.98e0 , units="(mol m-3)^-0.5 yr-1", + description="rate constant"), + ) + + stoich_redox_FeIII_H2S = PB.RateStoich( + PB.VarProp("redox_FeIII_H2S", "mol H2S yr-1", "sulphide consumption (+ve) by Fe+++ reduction", + attributes=(:calc_total=>true,)), + ((-1.0, "H2S"), (-8.0, "FeIIIOx"), (8.0, "FeII"), (1.0, "SO4"), (+14.0, "TAlk")), + sms_prefix="", + sms_suffix="_sms", + processname="redox", + ) + +end + +function PB.register_methods!(rj::ReactionRedoxFeIII_H2S) + + @info "register_methods! $(PB.fullname(rj))" + + vars = [ + rj.stoich_redox_FeIII_H2S.ratevartemplate, + PB.VarDep("H2S_conc", "mol m-3", "H2S concentration"), + PB.VarDep("FeIIIOx_conc", "mol m-3", "FeIIIOx concentration"), + PB.VarDep("volume", "m3", "box fluid volume"), + ] + + PB.add_method_do!( + rj, + do_FeIII_H2S_rate, + ( + PB.VarList_namedtuple(vars), + ); + ) + + PB.add_method_do!(rj, rj.stoich_redox_FeIII_H2S) + + PB.add_method_do_totals_default!(rj) + PB.add_method_initialize_zero_vars_default!(rj) + return nothing +end + + +function do_FeIII_H2S_rate(m::PB.ReactionMethod, pars, (vars, ), cellrange::PB.AbstractCellRange, deltat) + + rateparval = pars.R_FeIII_H2S[] + + for i in cellrange.indices + # mol yr-1 + vars.redox_FeIII_H2S[i] = ( rateparval # (mol m-3)^-0.5 yr-1 + * sqrt(max(PB.get_total(vars.H2S_conc[i]), 0.0)) # (mol m-3)^0.5 + * max(PB.get_total(vars.FeIIIOx_conc[i]), 0.0) # mol m-3 + * vars.volume[i]) # m3 + end + + return nothing +end + + +""" + ReactionRedoxFeS_O2 + +FeS oxidation + + FeS + 2 O2 -> Fe++ + SO4-- + +# Parameters +$(PARS) + +# Methods and Variables for default Parameters +$(METHODS_DO) +""" +Base.@kwdef mutable struct ReactionRedoxFeS_O2{P} <: PB.AbstractReaction + base::PB.ReactionBase + + pars::P = PB.ParametersTuple( + PB.ParDouble("R_FeS_O2", 1e5*1e-3 , units="(mol m-3)-1 yr-1", + description="rate constant"), + ) + + stoich_redox_FeS_O2 = PB.RateStoich( + PB.VarProp("redox_FeS_O2", "mol FeS yr-1", "FeS consumption (+ve) by oxidation by O2", + attributes=(:calc_total=>true,)), + ((-1.0, "FeS"), (-2.0, "O2"), (1.0, "FeII"), (1.0, "SO4"),), + sms_prefix="", + sms_suffix="_sms", + processname="redox", + ) + +end + +function PB.register_methods!(rj::ReactionRedoxFeS_O2) + + @info "register_methods! $(PB.fullname(rj))" + + ratevar = rj.stoich_redox_FeS_O2.ratevartemplate + concvars = [ + PB.VarDep("FeS_conc", "mol m-3", "FeS concentration"), + PB.VarDep("O2_conc", "mol m-3", "O2 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_FeS_O2), + name="do_redox_FeS_O2", + ) + + PB.add_method_do!(rj, rj.stoich_redox_FeS_O2) + + PB.add_method_do_totals_default!(rj) + PB.add_method_initialize_zero_vars_default!(rj) + return nothing +end + + +""" + ReactionRedoxFeS2pyr_O2 + +Pyrite oxidation + + FeS2 + 3.5 O2 + H2O -> Fe++ + 2 SO4-- + 2 H+ + +# Parameters +$(PARS) + +# Methods and Variables for default Parameters +$(METHODS_DO) +""" +Base.@kwdef mutable struct ReactionRedoxFeS2pyr_O2{P} <: PB.AbstractReaction + base::PB.ReactionBase + + pars::P = PB.ParametersTuple( + PB.ParDouble("R_FeS2pyr_O2", 1e3*1e-3 , units="(mol m-3)-1 yr-1", + description="rate constant"), + ) + + stoich_redox_FeS2pyr_O2 = PB.RateStoich( + PB.VarProp("redox_FeS2pyr_O2", "mol FeS2pyr yr-1", "pyrite consumption (+ve) by oxidation by O2", + attributes=(:calc_total=>true,)), + ((-1.0, "FeS2pyr"), (-3.5, "O2"), (1.0, "FeII"), (2.0, "SO4"), (-2.0, "TAlk")), + sms_prefix="", + sms_suffix="_sms", + processname="redox", + ) + +end + +function PB.register_methods!(rj::ReactionRedoxFeS2pyr_O2) + + @info "register_methods! $(PB.fullname(rj))" + + ratevar = rj.stoich_redox_FeS2pyr_O2.ratevartemplate + concvars = [ + PB.VarDep("FeS2pyr_conc", "mol m-3", "FeS2pyr (pyrite) concentration"), + PB.VarDep("O2_conc", "mol m-3", "O2 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_FeS2pyr_O2), + name="do_redox_FeS2pyr_O2", + ) + + PB.add_method_do!(rj, rj.stoich_redox_FeS2pyr_O2) + + PB.add_method_do_totals_default!(rj) + PB.add_method_initialize_zero_vars_default!(rj) + return nothing +end + + +###################################################################################################################### +# generic functions +###################################################################################################################### # 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)