diff --git a/docs/src/paleo_references.bib b/docs/src/paleo_references.bib index eabbf86..b1a1719 100644 --- a/docs/src/paleo_references.bib +++ b/docs/src/paleo_references.bib @@ -53,6 +53,35 @@ @article{Dale2015 year = {2015} } +@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{VanCappellen1996a, author = {{Van Cappellen}, Philippe and Wang, Y.}, doi = {10.2475/ajs.296.3.197}, diff --git a/examples/CONPSFeMn/PALEO_examples_sediment_NNFeMn.jl b/examples/CONPSFeMn/PALEO_examples_sediment_NNFeMn.jl new file mode 100644 index 0000000..cad5b12 --- /dev/null +++ b/examples/CONPSFeMn/PALEO_examples_sediment_NNFeMn.jl @@ -0,0 +1,124 @@ +using Logging +using Sundials +import LineSearches +using BenchmarkTools + +using Plots + +import PALEOboxes as PB +import PALEOmodel +import PALEOsediment + + +global_logger(ConsoleLogger(stderr,Logging.Info)) + +include("config_sediment_expts.jl") +include("../plot_sediment.jl") + +model = PB.create_model_from_config( + joinpath(@__DIR__, "PALEO_examples_sediment_NNFeMn_cfg.yaml"), + "sediment_Corg_O2NNMnFeS", +) + +############################# +# Steady state solutions +############################ + +config_sediment_expts(model, + [ + # oceanfloor oxygen + ("initial_value", "oceanfloor.O2_conc", [0.200, 0.001, 0.200, 0.001] ), + + # Corg input + # sensitivity test: remove most-reactive bin, maintaining same Corg input + # ("set_par", "sediment", "reservoir_Corg", "k_dist_modifier", [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0]), # remove 1 highest-k fractions + + # Fe & Mn input and remin + # high value 75 umol cm-2 yr-1 (0.75 mol m-2 yr-1, Van Cappellen 1996 test case) + # ("initial_value", "fluxOceanfloor.particulateflux_FeHR", 7.0.*[0.112, 0.112, 0.0177, 0.0177]), # mol m-2 yr-1 + # high value 40 um cm-2 yr-1 (0.4 mol m-2 yr-1, Van Cappellen 1996 test case) + # ("initial_value", "fluxOceanfloor.particulateflux_MnHR", 12.0.*[0.0342, 0.0342, 0.0101, 0.0101]), # mol m-2 yr-1 + # low values for limiting conc + # ("set_par", "sediment", "reminsed", "MnIVOxreminlimit", 40.0), # (mol m-3) default = 16e-6*2.5e6 (16 umol MnO2 g-1 assuming dry density is 2.5 g/cm^3 (= 2.5e6 g m^-3), Van Cappellen 1996) + # ("set_par", "sediment", "reminsed", "FeIIIOxreminlimit", 250.0), # (mol m-3) default = 100e-6*2.5e6 (100 umol Fe(OH)3 g-1 assuming dry density is 2.5 g/cm^3 (= 2.5e6 g m^-3), Van Cappellen 1996) + # ("set_par", "sediment", "reminsed", "FeIIIOxreminlimit", 0.1*250.0), # (mol m-3) = 100e-6*2.5e6 (100 umol Fe(OH)3 g-1 assuming dry density is 2.5 g/cm^3 (= 2.5e6 g m^-3), Van Cappellen 1996) + + + # bio rates + ("initial_value", "oceanfloor.alpha", [465.0, 465.0, 114.0, 114.0]), # yr-1 bioirrigation coefficient at surface + # ("initial_value", "oceanfloor.alpha", 0.5.*[465.0, 465.0, 114.0, 114.0]), # yr-1 reduced bioirrigation coefficient at surface + # ("initial_value", "oceanfloor.alpha", 0.0), # yr-1 no bioirrigation + + # Fe-S system + # Pyrite formation rate + # ("set_par", "sediment", "PyrH2S", "R_Pyr_H2S", 0.0), # zero rate (disable) pyrite formation + ("set_par", "sediment", "PyrH2S", "R_Pyr_H2S", 1e2), # 1e5 M yr-1 = 1e5*1e-3 (mol m-3) yr-1, Dale (2015) + # pyrite oxidation rate + ("set_par", "sediment", "redox_FeS2pyr_O2", "R_FeS2pyr_O2", 1.0), # (mol m-3)-1 yr-1, 1e3 M-1 yr-1 = 1e3*1e-3, Dale (2015) + # ("set_par", "sediment", "redox_FeS2pyr_O2", "R_FeS2pyr_O2", 100.0), # test x100 rate + # ("set_par", "sediment", "redox_FeS2pyr_O2", "R_FeS2pyr_O2", 0.0), # disable pyrite oxidation + ] +) + +tspan=(0.0, 100000.0) +# tspan=(0.0, 10.0) + +initial_state, modeldata = PALEOmodel.initialize!(model) + +run = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + +# PTC, Newton, no line search +# Bounds and max step size for Newton solve. NB requires min/max ratio for robustness so check all tracers initial_value is not zero +newton_min, newton_max, newton_min_ratio, newton_max_ratio = 1e-80, 1e6, 0.1, 10.0 +PALEOmodel.SteadyState.steadystate_ptcForwardDiff( + run, initial_state, modeldata, tspan, 1e-6, + deltat_fac=2.0, + solvekwargs=( + ftol=1e-7, + # ftol=1e-3, + # iterations=20, + iterations=50, # CFA formation has a discontinuous derivative -> need more iterations ? + method=:newton, + linesearch=LineSearches.Static(), + apply_step! = PALEOmodel.SolverFunctions.StepClampMultAll!(newton_min, newton_max, newton_min_ratio, newton_max_ratio), + # store_trace=true, + # show_trace=true, + ), + verbose=false, +) + + + +############################ +# Plot +############################ + +# single plots +# plotlyjs(size=(750, 565)) +# pager = PALEOmodel.DefaultPlotPager() + +# multiple plots per screen +colrange=1:4 +gr(size=(1200, 800)) +pager = PALEOmodel.PlotPager((1, 4), (legend_background_color=nothing, margin=(5, :mm))) + +plot_phi(run.output; colrange, pager) +plot_w(run.output; colrange, pager) +plot_biorates(run.output; colrange, pager) +plot_Corg_O2(run.output; Corgs=["Corg",], colT=[first(tspan), last(tspan)], colrange, pager) +plot_Corg_RCmultiG(run.output; colrange, pager) +plot_solutes(run.output; colT=[first(tspan), last(tspan)], solutes=["P", "DIC", "TAlk", "NO3", "NO2", "NH4", "SO4", "SmIIaqtot", "CH4", "H2", "MnII", "FeIIaqtot"], colrange, pager) +plot_FeP(run.output; colT=[first(tspan), last(tspan)], colrange, pager) +plot_sediment_FeS_summary(run.output; colrange, pager) +plot_solids(run.output; colT=[first(tspan), last(tspan)], solids=["MnHR", "MnMR", "FeHR", "FeMR", "FePR", "FeSm", "FeS2pyr", "PFeHR", "PFeMR", "CFA"], colrange, pager) +plot_carbchem(run.output; include_constraint_error=true, colT=last(tspan), colrange, pager) +plot_budget(run.output; name="P", solids=["Corg", "PFeHR", "PFeMR", "CFA"], stoich_factors=Dict("Corg"=>1/106), solutes=["P"], pager, colrange, concxscale=:log10, concxlims=(1e-3, 1e3)) +plot_budget(run.output; name="Mn", solids=["MnHR", "MnMR"], solutes=["MnII"], pager, colrange, concxscale=:log10, concxlims=(1e-3, 1e3)) +plot_budget(run.output; name="Fe", solids=["FeHR", "FeMR", "FePR", "FeSm", "FeS2pyr"], solutes=["FeIIaqtot"], pager, colrange, concxscale=:log10, concxlims=(1e-3, 1e4)) +plot_budget(run.output; name="S", solids=["FeSm", "FeS2pyr"], solutes=["SmIIaqtot", "SO4"], stoich_factors=Dict("FeS2pyr"=>2), pager, colrange, concxscale=:log10, concxlims=(1e-3, 1e4), fluxylims=(-1, 1) ) +plot_rates(run.output; colT=[first(tspan), last(tspan)], remin_rates=["reminOrgOxO2", "reminOrgOxNO2", "reminOrgOxNO3", "reminOrgOxMnIVOx", "reminOrgOxFeIIIOx", "reminOrgOxSO4", "reminOrgOxCH4"], colrange, pager) +plot_conc_summary(run.output; species=["O2", "NO3", "NO2", "NH4", "P", "SmIIaqtot", "FeIIaqtot", "MnII", "CH4"], pager, colrange, xscale=:log10, xlims=(1e-3, 1e0)) + +pager(:newpage) + + diff --git a/examples/CONPSFeMn/PALEO_examples_sediment_NNFeMn_cfg.yaml b/examples/CONPSFeMn/PALEO_examples_sediment_NNFeMn_cfg.yaml new file mode 100644 index 0000000..c32b72f --- /dev/null +++ b/examples/CONPSFeMn/PALEO_examples_sediment_NNFeMn_cfg.yaml @@ -0,0 +1,831 @@ +######################################################## +# Dale (2015) shelf/slope test cases with Corg, O2, NO2, NO3, Mn, Fe, SO4/H2S, P, with FeS system and pyrite +# Default config has four columns: 2*shelf, 2*slope +# Set num_columns and per-column parameters to define a different ensemble. +# NB: no N, Mn so results differ in detail from the paper. +############################################################## +sediment_Corg_O2NNMnFeS: + parameters: + num_columns: 4 + CIsotope: ScalarData + SIsotope: ScalarData + domains: + global: + # scalar domain + + fluxOceanfloor: + reactions: + floorstubparticulateflux: + class: ReactionFluxTarget + parameters: + const_stub: true + target_prefix: particulateflux_ + fluxlist: ["Corg", "MnHR", "MnMR", "FeHR", "FeMR", "FePR", "FeSm", "FeS2pyr", "PFeHR", "PFeMR", "CFA"] + + variable_attributes: + # Dale (2015) test cases + # shelf shelf slope slope + # + # 9.4 9.4 3 3 mmol/m^2/d + particulateflux_Corg:initial_value: [3.43, 3.43, 1.10, 1.10] # mol/column/yr-1 + # particulateflux_Corg:initial_delta: -25.0 + # 187 30 MnT umol/m^2/d, 1/2 of this to each of MnHR, MnMR + particulateflux_MnHR:initial_value: [0.0342, 0.0342, 0.0101, 0.0101] # mol yr-1 = umol/m^2/d * 1e-6 * area m^2 * 365.25 + particulateflux_MnMR:initial_value: [0.0342, 0.0342, 0.0101, 0.0101] # mol yr-1 = umol/m^2/d * 1e-6 * area m^2 * 365.25 + # 1840 290 FeT umol/m^2/d, 1/6 of this to each of FeHR, FeMR, FePR + particulateflux_FeHR:initial_value: [0.112, 0.112, 0.0177, 0.0177] # mol yr-1 = umol/m^2/d * 1e-6 * area m^2 * 365.25 + particulateflux_FeMR:initial_value: [0.112, 0.112, 0.0177, 0.0177] # mol yr-1 = umol/m^2/d * 1e-6 * area m^2 * 365.25 + particulateflux_FePR:initial_value: [0.112, 0.112, 0.0177, 0.0177] # mol yr-1 = umol/m^2/d * 1e-6 * area m^2 * 365.25 + particulateflux_FeSm:initial_value: 0.0 + particulateflux_FeS2pyr:initial_value: 0.0 + particulateflux_PFeHR:initial_value: 0.0 + particulateflux_PFeMR:initial_value: 0.0 + particulateflux_CFA:initial_value: 0.0 + + floorstubsoluteflux: + class: ReactionFluxTarget + parameters: + target_prefix: soluteflux_ + fluxlist: ["O2", "DIC", "TAlk", "P", "SO4", "SmIIaqtot", "CH4", "H2", "FeIIaqtot", "MnII", "NO2", "NO3", "NH4"] + + + fluxOceanBurial: + reactions: + + burialfluxes: + class: ReactionFluxTarget + parameters: + target_prefix: flux_ + fluxlist: ["Corg_1", "Corg_2", "Corg_3", "Corg_4", "Corg_5", "Corg_6", "Corg_7", "Corg_8", "Corg_9", "Corg_10", "Corg_11", "Corg_12", "Corg_13", "Corg_14", + "MnHR", "MnMR", "FeHR", "FeMR", "FePR", "FeSm", "FeS2pyr", "PFeHR", "PFeMR", "CFA"] + + total_Corg: # add up the multi-G fractions -> flux_Corg + class: ReactionVectorSum + parameters: + vars_prefix: flux_ + vars_to_add: ["Corg_1", "Corg_2", "Corg_3", "Corg_4", "Corg_5", "Corg_6", "Corg_7", "Corg_8", "Corg_9", "Corg_10", "Corg_11", "Corg_12", "Corg_13", "Corg_14" ] + variable_links: + sum: flux_Corg + + oceanfloor: + + reactions: + grid: + class: ReactionUnstructuredVectorGrid + parameters: + ncells: external%num_columns + + floorstubphys: + class: ReactionConst + parameters: + constnames: ["Afloor", "zfloor", "temp", "sal", "phi", "phimin", "w_accum", "zbioturb", "zbioirrig", "Dbio", "alpha"] + variable_attributes: + # Dale (2015) test cases + # shelf shelf slope slope + Afloor:initial_value: [1.0, 1.0, 1.0, 1.0] # m^2 column area + zfloor:initial_value: [-100.0, -100.0, -600.0, -600.0] # m + temp:initial_value: [283.15, 283.15, 281.15, 281.15] # K + sal:initial_value: 35.0 # psu + phi:initial_value: 0.9 # surface porosity + phimin:initial_value: 0.7 # minimum porosity + w_accum:initial_value: [100e-5, 100e-5, 16e-5, 16e-5] # m yr-1 (100, 16 cm kyr-1) + zbioturb:initial_value: 0.03 # m, 3 cm bioturbation characteristic depth + zbioirrig:initial_value: [0.02, 0.02, 0.0075, 0.0075] # m biorrigation characteristic depth + Dbio:initial_value: [28e-4, 28e-4, 18e-4, 18e-4] # m^2 yr-1 (28, 16 cm^2 yr-1) + alpha:initial_value: [465.0, 465.0, 114.0, 114.0] # yr-1 bioirrigation coefficient at surface + + floorstubsoluteconc: + class: ReactionConst + parameters: + constnames: ["O2_conc", "DIC_conc", "TAlk_conc", "P_conc", "SO4_conc", "SmIIaqtot_conc", "CH4_conc", "H2_conc", "FeIIaqtot_conc", "MnII_conc", "NO2_conc", "NO3_conc", "NH4_conc"] + variable_attributes: + # Dale (2015) test cases + # shelf shelf slope slope + O2_conc:initial_value: [0.200, 0.001, 0.200, 0.001] # (mol m-3) + DIC_conc:initial_value: 2.0 # (mol m-3) 2000 uM + TAlk_conc:initial_value: 2.2 # (mol m-3) 2200 uM + NO3_conc:initial_value: 35e-3 # 35 uM + NO2_conc:initial_value: 0.0 + NH4_conc:initial_value: 1e-3 # 1 uM + SO4_conc:initial_value: 28756.0e-3 # concentration mol m-3 ~ 28e-3 mol/kg * 1027 kg m-3 + SmIIaqtot_conc:initial_value: 0.0 # mol m-3 + CH4_conc:initial_value: 0.0 # mol m-3 + H2_conc:initial_value: 0.0 # mol m-3 + FeIIaqtot_conc:initial_value: 0.0 # mol m-3 + MnII_conc:initial_value: 0.0 + P_conc:initial_value: 0.001 # mol m-3 + + + transfer_particulatefluxOceanfloor: + class: ReactionFluxTransfer + parameters: + transfer_matrix: Identity + input_fluxes: fluxOceanfloor.particulateflux_$fluxname$ + output_fluxes: sediment.oceanfloor.$fluxname$_sms + variable_links: + + sediment: + reactions: + bioratesed: + class: ReactionSedimentBioRates + parameters: + f_bioTurbDepth: Exp2Cutoff + f_bioIrrigDepth: Exp1Cutoff + separate_zbio: true + f_bioO2: Dale2015 + bioO2halfmax: 20e-3 # (mol m-3) (Dale 2015 'a', 20 uM) + bioO2decreaserate: 12e-3 # (mol m-3) (Dale 2015 'b', 12 uM) + + transportsed: + class: ReactionSedimentTransport + parameters: + L: 0.30 # (m) column length + ncellspercol: 100 + f_grid: quadratic + grid_eta: 0.075 + zdbl: 0.04e-2 # diffusive boundary layer at sediment-water interface + f_porosity: ExpAtten # + zpor: 0.1 # (m) lengthscale for porosity decrease + w_solute: true # include solute advection (w_solute = w_solid at base of column) + + reservoir_O2: + class: ReactionReservoirTotal + variable_links: + R*: O2* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 200e-3 + R_conc:diffusivity_speciesname: O2 + + reservoir_NO2: + class: ReactionReservoirTotal + variable_links: + R*: NO2* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 1e-3 + R_conc:diffusivity_speciesname: NO2 + + reservoir_NO3: + class: ReactionReservoirTotal + variable_links: + R*: NO3* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 10e-3 + R_conc:diffusivity_speciesname: NO3 + + reservoir_NH4: + class: ReactionReservoirTotal + variable_links: + R*: NH4* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 1e-3 + R_conc:diffusivity_speciesname: NH4 + + reservoir_SO4: + class: ReactionReservoirTotal + variable_links: + R*: SO4* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 20000e-3 + R_conc:diffusivity_speciesname: SO4 + + reservoir_SmIIaqtot: + class: ReactionReservoirTotal + variable_links: + R*: SmIIaqtot* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 200e-3 + R_conc:diffusivity_speciesname: H2S + + reservoir_CH4: + class: ReactionReservoirTotal + variable_links: + R*: CH4* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 200e-3 + R_conc:diffusivity_speciesname: CH4 + + reservoir_H2: + class: ReactionReservoirTotal + parameters: + + variable_links: + R*: H2* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration mol m-3 ~ 1e-9 mol/kg * 1027 kg m-3 + # R:initial_delta: 0.0 + R:norm_value: 1.0e-3 # for scaling only + R_conc:diffusivity_speciesname: H2 + + reservoir_P: + class: ReactionReservoirTotal + variable_links: + R*: P* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 2e-3 + R_conc:diffusivity_speciesname: PO4 + # R_conc:gamma: 0.2 # reduced bioirrigation + + reservoir_DIC: + class: ReactionReservoirTotal + parameters: + field_data: external%CIsotope + variable_links: + R*: DIC* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 1000e-3 + R_conc:diffusivity_speciesname: HCO3 # use for molecular diffusivity HCO3- as the major species + + reservoir_TAlk: + class: ReactionReservoirTotal + variable_links: + R*: TAlk* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 1000e-3 + R_conc:diffusivity_speciesname: Na # TODO - use major cations for molecular diffusivity + + reservoir_FeIIaqtot: + class: ReactionReservoirTotal + variable_links: + R*: FeIIaqtot* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 2e-3 + R_conc:diffusivity_speciesname: Fe + R_conc:gamma: 0.2 # reduced bioirrigation for FeII + + reservoir_MnII: + class: ReactionReservoirTotal + variable_links: + R*: MnII* + volume: volume_solute + variable_attributes: + R_conc:vphase: VP_Solute + R:initial_value: 1e-12 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw) + R:norm_value: 1e-3 + R_conc:diffusivity_speciesname: Mn + + reservoir_Corg: + class: ReactionRCmultiG + parameters: + # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 + k_bin_edges: [1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0, 1e1, 1e2] # yr-1 + k_dist_modifier: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] + a: 3e-4 # yr-1 reactive-continuum model average lifetime + v: 0.125 # reactive-continuum model shape parameter + field_data: external%CIsotope + + variable_links: + volume: volume_solid + # POC_decay: DIC_sms # dump all Corg -> DIC for test case + + variable_attributes: + Corg_*_conc:vphase: VP_Solid + Corg_*:initial_value: 1e-12 + + reservoir_MnHR: + class: ReactionReservoirTotal + variable_links: + R*: MnHR* + volume: volume_solid + variable_attributes: + R_conc:vphase: VP_Solid + R:initial_value: 1e-12 # concentration m-3 solid phase + # mol m-3 = g cm-3 / (g mol-1) * m^3 cm-3 + R:norm_value: 1.0 # mol m-3 solid + + reservoir_MnMR: + class: ReactionReservoirTotal + variable_links: + R*: MnMR* + volume: volume_solid + variable_attributes: + R_conc:vphase: VP_Solid + R:initial_value: 1e-12 # concentration m-3 solid phase + # mol m-3 = g cm-3 / (g mol-1) * m^3 cm-3 + R:norm_value: 1.0 # mol m-3 solid + + reservoir_FeHR: + class: ReactionReservoirTotal + variable_links: + R*: FeHR* + volume: volume_solid + variable_attributes: + R_conc:vphase: VP_Solid + R:initial_value: 1e-12 # concentration m-3 solid phase + # mol m-3 = g cm-3 / (g mol-1) * m^3 cm-3 + R:norm_value: 1.0 # mol m-3 solid + + reservoir_FeMR: + class: ReactionReservoirTotal + variable_links: + R*: FeMR* + volume: volume_solid + variable_attributes: + R_conc:vphase: VP_Solid + R:initial_value: 1e-12 # concentration m-3 solid phase + # mol m-3 = g cm-3 / (g mol-1) * m^3 cm-3 + R:norm_value: 1.0 # mol m-3 solid + + reservoir_FePR: + class: ReactionReservoirTotal + variable_links: + R*: FePR* + volume: volume_solid + variable_attributes: + R_conc:vphase: VP_Solid + R:initial_value: 1e-12 # concentration m-3 solid phase + # mol m-3 = g cm-3 / (g mol-1) * m^3 cm-3 + R:norm_value: 1.0 # mol m-3 solid + + reservoir_FeSm: + class: ReactionReservoirTotal + + variable_links: + R*: FeSm* + volume: volume_solid + variable_attributes: + R_conc:vphase: VP_Solid + R:initial_value: 1e-12 # concentration mol m-3 ~ 1e-9 mol/kg * 1027 kg m-3 + R:norm_value: 1.0e-3 # for scaling only + + reservoir_FeS2pyr: + class: ReactionReservoirTotal + + variable_links: + R*: FeS2pyr* + volume: volume_solid + variable_attributes: + R_conc:vphase: VP_Solid + R:initial_value: 1e-12 # concentration mol m-3 ~ 1e-9 mol/kg * 1027 kg m-3 + R:norm_value: 1.0e-3 # for scaling only + + reservoir_PFeHR: # FeHR-bound phosphorus + class: ReactionReservoirTotal + variable_links: + R*: PFeHR* + volume: volume_solid + variable_attributes: + R_conc:vphase: VP_Solid + R:initial_value: 1e-12 # concentration m-3 solid phase + # mol m-3 = g cm-3 / (g mol-1) * m^3 cm-3 + R:norm_value: 1.0 # mol m-3 solid + + reservoir_PFeMR: # FeMR-bound phosphorus + class: ReactionReservoirTotal + variable_links: + R*: PFeMR* + volume: volume_solid + variable_attributes: + R_conc:vphase: VP_Solid + R:initial_value: 1e-12 # concentration m-3 solid phase + # mol m-3 = g cm-3 / (g mol-1) * m^3 cm-3 + R:norm_value: 1.0 # mol m-3 solid + + reservoir_CFA: # carbonate fluorapatite + class: ReactionReservoirTotal + variable_links: + R*: CFA* + volume: volume_solid + variable_attributes: + R_conc:vphase: VP_Solid + R:initial_value: 1e-12 # concentration m-3 solid phase + # mol m-3 = g cm-3 / (g mol-1) * m^3 cm-3 + R:norm_value: 1.0 # mol m-3 solid + + reservoir_Ca: + class: ReactionReservoirConst + variable_links: + R*: Ca* + variable_attributes: + R_conc:initial_value: 10.56 # mol m-3 10.56 mM modern Ca concentration + + reservoir_B: + class: ReactionReservoirConst + variable_links: + R*: B* + variable_attributes: + R_conc:initial_value: 0.427 # mol m-3 0.427 mM modern B concentration + + reservoir_F: + class: ReactionReservoirConst + variable_links: + R*: F* + variable_attributes: + R_conc:initial_value: 70.17e-3 # mol m-3 70.17 uM modern F concentration + + carbchem: + class: ReactionCO2SYS + parameters: + # TODO NH4 ? + components: ["Ci", "B", "S", "F", "Omega"] # , "H2S"] no H2S as that is handled by ReactionFeSaq + defaultconcs: [] # ["TF", "TB", "Ca"] + solve_pH: constraint # H as a state variable + outputs: ["H", "CO2", "HCO3", "CO3", "OmegaCA", "OmegaAR"] + + variable_links: + volume: volume_solute + H: H_conc + TCi_conc: DIC_conc + TS_conc: SO4_conc + TAlk_calc: TAlk_calc + # TH2S_conc: H2S_conc + TB_conc: B_conc + TF_conc: F_conc + CO2: CO2_conc + HCO3: HCO3_conc + CO3: CO3_conc + OmegaCA: OmegaCA + OmegaAR: OmegaAR + + # reservoir_Hp: + # # for testing: fixed H+ for constant pH 8 + # class: ReactionReservoirConst + # variable_links: + # R*: Hp* + # variable_attributes: + # R_conc:initial_value: 1e-5 # mol m-3 10^-8.0 *1e3 [H+] at pH 8.0 + + FeSaq: + class: ReactionFeSaq + parameters: + add_TAlk_calc: true # add S-II and FeII contributions to TAlk_calc + variable_links: + volume: volume_solute + Hp_conc: H_conc # defined by carbonate chemistry solver + + + pocdecaycomponents: + class: ReactionFluxToComponents + parameters: + outputflux_prefix: remin_ + outputflux_names: ["Corg", "N", "P"] + # 106/106 16/106 1/106 + outputflux_stoich: [1.0, 0.15094, 0.009434] # must match input Corg stoich + variable_links: + inputflux: POC_decay + + MnHRdecay: + class: ReactionParticleDecay + parameters: + decay_timescale: 1.7 # yr, Dale (2015) + variable_links: + Particle*: MnHR* + decayflux: MnMR_sms + + FeHRdecay: + class: ReactionParticleDecay + parameters: + decay_timescale: 0.7 # yr, Dale (2015) + show_decayrate: true # needed for FeP coprecipitation + variable_links: + Particle*: FeHR* + decayflux: FeMR_sms + decayrate: rate_FeHR_FeMR + + reminsed: + class: ReactionReminO2_NN_Mn_Fe_SO4_CH4 + + parameters: + oxreminlimit: 1e-3 # (mol m-3) Dale (2015) + NO2reminlimit: 10e-3 # (mol m-3) Dale (2015) + NO3reminlimit: 10e-3 # (mol m-3) Dale (2015) + MnIVOxreminlimit: 40.0 # (mol m-3) = 16e-6*2.5e6 (16 umol MnO2 g-1 assuming dry density is 2.5 g/cm^3 (= 2.5e6 g m^-3), Van Cappellen 1996) + FeIIIOxreminlimit: 250.0 # (mol m-3) = 100e-6*2.5e6 (100 umol Fe(OH)3 g-1 assuming dry density is 2.5 g/cm^3 (= 2.5e6 g m^-3), Van Cappellen 1996) + SO4reminlimit: 500.0e-3 # (mol m-3) Dale (2015) + + variable_links: + FeIIIOx_conc: FeHR_conc + MnIVOx_conc: MnHR_conc + soluteflux_*: "*_sms" + soluteflux_H2S: SmIIaqtot_sms + soluteflux_FeII: FeIIaqtot_sms + soluteflux_FeIIIOx: FeHR_sms + soluteflux_MnIVOx: MnHR_sms + soluteflux_TNH3: NH4_sms + + redox_NH4_NO2: + class: ReactionRedoxNH4_NO2 + parameters: + R_NH4_NO2: 1e5 # (mol m-3)-1 1e8*1e-3 (Dale 2015) + variable_links: + volume: volume_solute + + redox_NH4_O2: + class: ReactionRedoxNH4_O2 + parameters: + R_NH4_O2: 1e4 # (mol m-3)-1 1e7*1e-3 (Dale 2015) + variable_links: + volume: volume_solute + + redox_NO2_O2: + class: ReactionRedoxNO2_O2 + parameters: + R_NO2_O2: 1e4 # (mol m-3)-1 1e7*1e-3 (Dale 2015) + variable_links: + volume: volume_solute + + redox_H2S_O2: + class: ReactionRedoxH2S_O2 + parameters: + R_H2S_O2: 3e5 # Boudreau (1997) + variable_links: + volume: volume_solute + H2S_conc: H2Stot_conc + H2S_sms: SmIIaqtot_sms + + redox_CH4_O2: + class: ReactionRedoxCH4_O2 + + parameters: + # R_CH4_O2: 1e4 # reduce rate for numerical stability + R_CH4_O2: 1e7 # (mol m-3) yr-1 (1e10 mol/l yr-1) Boudreau (1997) + variable_links: + volume: volume_solute + + redox_CH4_SO4: + class: ReactionRedoxCH4_SO4 + + parameters: + R_CH4_SO4: 10.0 # (mol m-3) yr-1 + + variable_links: + volume: volume_solute + H2S_sms: SmIIaqtot_sms + + redox_H2_O2: + class: ReactionRedoxH2_O2 + + parameters: + R_H2_O2: 1e6 # (mol m-3)-1 yr-1 (1e6 M-1 yr-1 = 1e6*1e-3 (mol m-3)-1 yr-1, Dale (2015). NB: arbitrary, just need to scavenge H2 produced during pyrite formation) + variable_links: + volume: volume_solute + + redox_H2_SO4: + class: ReactionRedoxH2_SO4 + + parameters: + R_H2_SO4: 1e6 # (mol m-3)-1 yr-1 (1e6 M-1 yr-1 = 1e6*1e-3 (mol m-3)-1 yr-1, Dale (2015). NB: arbitrary, just need to scavenge H2 produced during pyrite formation) + + variable_links: + volume: volume_solute + H2S_sms: SmIIaqtot_sms + + redox_MnII_O2: + class: ReactionRedoxMnII_O2 + parameters: + R_MnII_O2: 5e3 # (mol m-3)-1 yr-1, 5e6 M-1 yr-1 Dale (2015) + variable_links: + volume: volume_solute + MnIVOx*: MnHR* + + redox_FeII_O2: + class: ReactionRedoxFeII_O2 + parameters: + R_FeII_O2: 5e5 # (mol m-3)-1 yr-1, 5e8 M-1 yr-1 = 5e8*1e-3 (mol m-3)-1 yr-1, Dale (2015) + variable_links: + volume: volume_solute + FeII_sms: FeIIaqtot_sms + FeIIIOx*: FeHR* + + redox_FeII_NO3: + class: ReactionRedoxFeII_NO3 + parameters: + R_FeII_NO3: 1e2 # (mol m-3)-1 yr-1, 1e5 M-1 yr-1 = 1e5*1e-3 (mol m-3)-1 yr-1, Dale (2015) + variable_links: + volume: volume_solute + FeII_sms: FeIIaqtot_sms + FeIIIOx*: FeHR* + + redox_MnHR_FeII: + class: ReactionRedoxMnIV_FeII + parameters: + R_MnIV_FeII: 1e4 # (mol m-3 yr-1), Dale (2015) 1e7, 1e5 M-1 yr-1 for MnHR, MnMR (-> 1e4, 1e2 (mol m-3)-1 yr-1) + variable_links: + volume: volume_solid + redox_MnIV_FeII*: redox_MnHR_FeII* + MnIVOx*: MnHR* + FeII_sms: FeIIaqtot_sms + FeIIIOx*: FeHR* + + redox_MnMR_FeII: + class: ReactionRedoxMnIV_FeII + parameters: + R_MnIV_FeII: 1e2 # (mol m-3 yr-1), Dale (2015) 1e7, 1e5 M-1 yr-1 for MnHR, MnMR (-> 1e4, 1e2 (mol m-3)-1 yr-1) + variable_links: + volume: volume_solid + redox_MnIV_FeII*: redox_MnMR_FeII* + MnIVOx*: MnMR* + FeII_sms: FeIIaqtot_sms + FeIIIOx*: FeHR* + + redox_MnHR_H2S: + class: ReactionRedoxMnIV_H2S + + parameters: + R_MnIV_H2S: 1e2 # (mol m-3)^1 yr-1, Dale (2015) 1e5, 1e3 M-1 yr-1 for MnHR, MnMR (-> 1e2, 1e0 (mol m-3)-1 yr-1) + + variable_links: + volume: volume_solid + redox_MnIV_H2S*: redox_MnHR_H2S* + MnIVOx*: MnHR* + H2S_conc: H2Stot_conc + H2S_sms: SmIIaqtot_sms + + redox_MnMR_H2S: + class: ReactionRedoxMnIV_H2S + + parameters: + R_MnIV_H2S: 1e0 # (mol m-3)^1 yr-1, Dale (2015) 1e5, 1e3 M-1 yr-1 for MnHR, MnMR (-> 1e2, 1e0 (mol m-3)-1 yr-1) + + variable_links: + volume: volume_solid + redox_MnIV_H2S*: redox_MnMR_H2S* + MnIVOx*: MnMR* + H2S_conc: H2Stot_conc + H2S_sms: SmIIaqtot_sms + + redox_FeHR_H2S: + class: ReactionRedoxFeIII_H2S + + parameters: + R_FeIII_H2S: 3.16 # (mol m-3)^-0.5 yr-1, 100 M-0.5 yr-1 = 100*3.16e-2 (mol m-3)^-0.5 yr-1 + + variable_links: + volume: volume_solid + redox_FeIII_H2S*: redox_FeHR_H2S* + FeIIIOx*: FeHR* + FeII_sms: FeIIaqtot_sms + H2S_conc: H2Stot_conc + H2S_sms: SmIIaqtot_sms + + redox_FeMR_H2S: + class: ReactionRedoxFeIII_H2S + + parameters: + R_FeIII_H2S: 3.16e-3 # (mol m-3)^-0.5 yr-1, 0.1 M-0.5 yr-1 = 0.1*3.16e-2 (mol m-3)^-0.5 yr-1 + + variable_links: + volume: volume_solid + redox_FeIII_H2S*: redox_FeMR_H2S* + FeIIIOx*: FeMR* + FeII_sms: FeIIaqtot_sms + H2S_conc: H2Stot_conc + H2S_sms: SmIIaqtot_sms + + redox_FePR_H2S: + class: ReactionRedoxFeIII_H2S + + parameters: + R_FeIII_H2S: 9.49e-6 # (mol m-3)^-0.5 yr-1, 3e-4 M-0.5 yr-1 = 3e-4*3.16e-2 (mol m-3)^-0.5 yr-1 + + variable_links: + volume: volume_solid + redox_FeIII_H2S*: redox_FePR_H2S* + FeIIIOx*: FePR* + FeII_sms: FeIIaqtot_sms + H2S_conc: H2Stot_conc + H2S_sms: SmIIaqtot_sms + + redox_FeS_O2: + class: ReactionRedoxFeS_O2 + + parameters: + R_FeS_O2: 1e2 # (mol m-3)-1 yr-1, 1e5 M-1 yr-1 = 1e5*1e-3 + + variable_links: + volume: volume_solid + FeS_*: FeSm_* + FeII_sms: FeIIaqtot_sms + + redox_FeS2pyr_O2: + class: ReactionRedoxFeS2pyr_O2 + + parameters: + R_FeS2pyr_O2: 1.0 # (mol m-3)-1 yr-1, 1e3 M-1 yr-1 = 1e3*1e-3 + + variable_links: + volume: volume_solid + FeII_sms: FeIIaqtot_sms + + FeSm: + class: ReactionFeSm + + parameters: + rate_FeS_prec: 40.0 # mol m-3 yr-1 rate of FeS precipitation (van Cappellen & Wang value 1.5e-5 mol g yr-1 ~ 1.5-5 * 2.5e6 ~ 37.5 for rho = 2.5e6 g m-3) + rate_FeS_diss: 1e3 # yr-1 rate of FeS dissolution + + PyrH2S: + class: ReactionPyrH2S + + parameters: + R_Pyr_H2S: 1e2 # (mol m-3) yr-1, 1e5 M yr-1 = 1e5*1e-3 (mol m-3) yr-1, Dale (2015) + variable_links: + volume: volume_solid + H2Ssp_conc: H2Stot_conc # Dale (2015) write this for total [H2S] = [H2S] + [HS-], not H2S species + + PFeHR_coprecip: + class: ReactionPACoPrecip + + parameters: + A_rate_stoich_factors: [1.0, 1.0, 2.0, 2.0] + gamma: 0.1 + P_limit: 10e-3 # mol m-3 (10 uM, Dale 2016) + + variable_links: + rate_PA_coprecip*: rate_PFeHR_coprecip* + A_formation_rate_1: redox_FeII_O2 # mol Fe yr-1 + A_formation_rate_2: redox_FeII_NO3 # mol Fe yr-1 + A_formation_rate_3: redox_MnHR_FeII # mol Mn yr-1 + A_formation_rate_4: redox_MnMR_FeII # mol Mn yr-1 + PA_sms: PFeHR_sms + + PFeHR_release: + class: ReactionPARelease + + parameters: + A_rate_stoich_factors: [-4.0, 1.0] + + variable_links: + A_conc: FeHR_conc + PA_conc: PFeHR_conc + PA_theta: PFeHR_theta + rate_PA_release*: rate_PFeHR_release* + A_destruction_rate_1: reminOrgOxFeIIIOx # mol O2eq yr-1 + A_destruction_rate_2: redox_FeHR_H2S # mol Fe yr-1 + PA_sms: PFeHR_sms + + PFeHR_PFeMR: # transfer adsorbed P from FeHR to FeMR phase + class: ReactionPARelease + + parameters: + A_rate_stoich_factors: [1.0] + partition_frac: 0.25 # 1/4 of release P from FeHR -> FeMR, 3/4 to dissolved P + + variable_links: + A_conc: FeHR_conc + PA_conc: PFeHR_conc + PA_theta: PFeHR_theta_unused # not needed + rate_PA_release*: rate_PFeHR_PFeMR* + A_destruction_rate_1: rate_FeHR_FeMR + PA_sms: PFeHR_sms + P2_sms: PFeMR_sms # transfer partition_frac of released P to FeMR + + PFeMR_release: + class: ReactionPARelease + + parameters: + A_rate_stoich_factors: [1.0] + + variable_links: + A_conc: FeMR_conc + PA_conc: PFeMR_conc + PA_theta: PFeMR_theta + rate_PA_release*: rate_PFeMR_release* + A_destruction_rate_1: redox_FeMR_H2S # mol Fe yr-1 + PA_sms: PFeMR_sms + + CFA_formation: + class: ReactionCFAsimple + + parameters: + R_CFA: 1.0 # yr-1 CFA formation rate constant Dale etal (2016) quote 1.0 yr-1 + Peq_conc: 10e-3 # mol m-3 Dale etal (2016) quote 10e-6*1e3 + + sedimentfloor: + reactions: + + diff --git a/examples/CONPSFeMn/README.md b/examples/CONPSFeMn/README.md new file mode 100644 index 0000000..1b6327a --- /dev/null +++ b/examples/CONPSFeMn/README.md @@ -0,0 +1,35 @@ +# Sediment 1D reaction-transport configuration with C, O, N, P, S, Fe, Mn + + julia> include("PALEO_examples_sediment_NNFeMn.jl") + +Configuration is similar to shelf and slope test cases in +Dale (2015) GBC , with: +- a discrete (multi-G) approximation to POC reactive-continuum distribution +- NO2 and NO3 oxidants +- Mn oxides with two reactivity fractions +- Fe oxides with three reactivity fractions + +There are four sediment columns: oxic and anoxic 'shelf', oxic and anoxic 'slope' + +Differences from Dale (2015): + +- Fe-S system: + - solves speciation between FeII, H2S and FeSaq + (Rickard & Luther 2007 ) + with FeSm preciptiation based on FeSaq solubility threshold. + - Pyrite formation only includes 'Berzelius' pathway + (FeSm + H2S -> FeS2pyr + H2) + - No S0 + +- Includes a minimal model for P diagenesis including: + - P coprecipitation on FeHR, FeMR iron oxides + - CFA formation + + See: + - Slomp (1996) + - Reed (2011) + - Dale (2016) + - Zhao (2020) + + + diff --git a/examples/CONPSFeMn/config_sediment_expts.jl b/examples/CONPSFeMn/config_sediment_expts.jl new file mode 100644 index 0000000..2aabdfe --- /dev/null +++ b/examples/CONPSFeMn/config_sediment_expts.jl @@ -0,0 +1,38 @@ + +"test cases and examples for sediment" +function config_sediment_expts( + model, expts; +) + + for expt in expts + println("Add expt: ", expt) + if expt == "baseline" + # defaults + + elseif length(expt) == 5 && expt[1] == "set_par" + # generic parameter set (set_par, , , , rate_var / volume_var + +Convert per-cell rate (mol yr-1) to `rate_var / volume_var` (mol m-3 yr-1) +""" +function normalize_by_volume(rate_var::PALEOmodel.FieldArray, volume_var::PALEOmodel.FieldArray) + return PALEOmodel.FieldArray(rate_var.name*" / volume", rate_var.values./volume_var.values, rate_var.dims, rate_var.attributes) +end + function plot_rates( output; remin_rates=["reminOrgOxO2", "reminOrgOxSO4", "reminOrgOxCH4"], - colT=[-1e12, 0.1, 1.0, 10.0, 100.0, 1000.0, 1e12], + colT=1e12, colrange=1:3, pager=PALEOmodel.DefaultPlotPager(), ) for icol in colrange - pager( - plot( - title="Corg remin rate $icol", - output, - ["sediment.remin_Corg"], - (tmodel=colT, column=icol), - xlabel="remin (mol Corg yr-1)", - swap_xy=true - ) - ) + p = plot() + # remin_Corg at each time + for t in colT + volume = PALEOmodel.get_array(output, "sediment.volume", (tmodel=t, column=icol)) + remin_Corg = PALEOmodel.get_array(output, "sediment.remin_Corg", (tmodel=t, column=icol)) + remin_Corg_norm = normalize_by_volume(remin_Corg, volume) + plot!(p, remin_Corg_norm; swap_xy=true) + end + + plot!(p, title="Corg remin rate $icol", xlabel="remin (mol Corg m-3 yr-1)") + + pager(p) end for icol in colrange - pager( - plot( - title="Org matter ox rate $icol", - output, - "sediment.".*remin_rates, - (tmodel=colT[end], column=icol), # only last time - xlabel="ox rate (mol O2eq yr-1)", - swap_xy=true, - ) - ) + # remin_rates at last time + volume = PALEOmodel.get_array(output, "sediment.volume", (tmodel=last(colT), column=icol)) + p = plot() + + for r in remin_rates + rate = PALEOmodel.get_array(output, "sediment.".*r, (tmodel=last(colT), column=icol)) + rate_norm = normalize_by_volume(rate, volume) + plot!(p, rate_norm; swap_xy=true) + end + plot!(p; title="Org matter ox rate $icol", xlabel="ox rate (mol O2eq m-3 yr-1)") + + pager(p) end end @@ -307,4 +345,103 @@ function plot_carbchem( pager(:newpage) return nothing -end \ No newline at end of file +end + +function plot_FeP( + output; + FeP=["PFeHR_theta", "PFeMR_theta"], + colT=[-1e12, 0.1, 1.0, 10.0, 100.0, 1000.0, 1e12], + colrange=1:3, + pager=PALEOmodel.DefaultPlotPager(), +) + + for fp in FeP + for icol in colrange + pager(plot(title="[$fp] $icol", output, ["sediment.$(fp)"], (tmodel=colT, column=icol), xlabel="$fp P:Fe (mol/mol)", swap_xy=true)) + end + end + + pager(:newpage) + + return nothing +end + + +function plot_budget( + output; + name="Mn", + solids=["MnHR", "MnMR"], + solutes=["MnII"], + stoich_factors=Dict(), # eg Dict("Corg"=>1/106) for P:Corg ratio + concxscale=:log10, + concxlims=(1e-3, Inf), + fluxyscale=:identity, + fluxylims=(-Inf, Inf), + colT=last(PALEOmodel.get_array(output, "global.tmodel").values), # model time for column plots + colrange=1:3, + pager=PALEOmodel.DefaultPlotPager(), +) + # scale by stoich_factors[species] + function scale_stoich(Cconc, species, label) + if haskey(stoich_factors, species) + sf = stoich_factors[species] + Cconc *= sf + label *= " * $(round(sf, sigdigits=3))" + end + return Cconc, label + end + + # plot concentrations + all_species=vcat(solutes, solids) + for icol in colrange + p = plot(title="$name concentrations $icol") + for sp in all_species + conc = PALEOmodel.get_array(output, "sediment.$(sp)_conc", (column=icol, tmodel=colT)) + conc, label = scale_stoich(conc, sp, sp) + plot!(p, conc; label, swap_xy=true) + end + plot!(p, xlabel="conc (mol m-3)", xscale=concxscale, xlims=concxlims) + pager(p) + end + + # plot input/output fluxes and total budget + for icol in colrange + tmodel = PALEOmodel.get_array(output, "sediment.tmodel").values + total_flux_into_sediment = zeros(length(tmodel)) + p = plot(title="$name budget $icol") + for s in solutes + # +ve is sediment -> ocean + fluxname = "fluxOceanfloor.soluteflux_$s" + solute_flux = PALEOmodel.get_array(output, fluxname, (cell=icol,)) + solute_flux, label = scale_stoich(solute_flux, s, "-"*fluxname) + total_flux_into_sediment .-= solute_flux.values + plot!(p, -1*solute_flux; label, linestyle=:dash) + end + for s in solids + # +ve is into sediment + fluxname = "fluxOceanfloor.particulateflux_$s" + particulate_flux = PALEOmodel.get_array(output, fluxname, (cell=icol,)) + particulate_flux, label = scale_stoich(particulate_flux, s, fluxname) + total_flux_into_sediment .+= particulate_flux.values + plot!(p, particulate_flux; label) + end + for s in solids + # +ve is out of sediment + fluxname = "fluxOceanBurial.flux_$s" + burial_flux = PALEOmodel.get_array(output, fluxname, (cell=icol,)) + burial_flux, label = scale_stoich(burial_flux, s, "-"*fluxname) + total_flux_into_sediment .-= burial_flux.values + plot!(p, -1*burial_flux; label) + end + + plot!(p, tmodel, total_flux_into_sediment, label="total", color=:black) + + plot!(p, ylabel="$name flux (mol yr-1)", yscale=fluxyscale, ylims=fluxylims) + pager(p) + end + + pager(:newpage) + + return nothing +end +