Skip to content

Commit

Permalink
Merge pull request #2 from PALEOtoolkit/S_isotopes_example
Browse files Browse the repository at this point in the history 
Add  sulphur isotope example
  • Loading branch information
@sjdaines
sjdaines authored Feb 29, 2024
2 parents ca5dabe + eca8644 commit 87d34fb
Show file tree
Hide file tree
Showing 5 changed files with 133 additions and 12 deletions.
1 change: 1 addition & 0 deletions examples/Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
[deps]
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
Expand Down
83 changes: 83 additions & 0 deletions examples/boudreau1996/PALEO_examples_sediment_Sisotopes.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
using Logging
using Sundials
import LineSearches
using BenchmarkTools

using Plots

import PALEOboxes as PB
import PALEOmodel
import PALEOaqchem
import PALEOsediment

do_benchmarks = false

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_cfg.yaml"),
"sediment_Corg_O2";
modelpars=Dict("SIsotope"=>"IsotopeLinear"), # enable d34S isotopes
)

#############################
# Steady state solutions
############################

config_sediment_expts(model,
[
("initial_value", "oceanfloor.SO4_conc", 1000e-3), # 1mM SO4
]
)

tspan=(0.0, 10000.0)

initial_state, modeldata = PALEOmodel.initialize!(model)

run = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory())

# PTC, Newton, no line search
# TODO This form of Newton regularization doesn't work with isotopes !! (as isotope mol*delta state variables can be -ve)
# Bounds and max step size for Newton solve. NB some tracers start at zero so set newton_max_ratio=Inf
# newton_min, newton_max, newton_min_ratio, newton_max_ratio = 1e-80, Inf, 0.1, Inf
PALEOmodel.SteadyState.steadystate_ptcForwardDiff(
run, initial_state, modeldata, tspan, 1e-3,
deltat_fac=2.0,
solvekwargs=(
ftol=1e-7,
iterations=20,
method=:newton,
linesearch=LineSearches.Static(),
# TODO Newton regularization doesn't work with isotopes
# 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
gr(size=(900, 600))
pager = PALEOmodel.PlotPager((1,3), (legend_background_color=nothing, ))

plot_Corg_O2(run.output; Corgs=["Corg1", "Corg2"], colT=[first(tspan), last(tspan)], pager=pager)
plot_solutes(run.output; colT=[first(tspan), last(tspan)], pager=pager)
plot_solute_deltas(run.output; solutes=["SO4", "H2S"], colT=[first(tspan), last(tspan)], pager=pager)
plot_rates(run.output; colT=[first(tspan), last(tspan)], pager=pager)

pager(:newpage)


28 changes: 17 additions & 11 deletions examples/boudreau1996/PALEO_examples_sediment_cfg.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ sediment_Corg_O2:
class: ReactionFluxTarget
parameters:
target_prefix: soluteflux_
fluxlist: ["O2", "P", "SO4", "H2S", "CH4"]
fluxlist: ["O2", "P", "SO4::SIsotope", "H2S::SIsotope", "CH4"]


fluxOceanBurial:
Expand Down Expand Up @@ -79,13 +79,15 @@ sediment_Corg_O2:
floorstubsoluteconc:
class: ReactionConst
parameters:
constnames: ["O2_conc", "P_conc", "SO4_conc", "H2S_conc", "CH4_conc"]
constnames: ["O2_conc", "P_conc", "SO4_conc::SIsotope", "H2S_conc::SIsotope", "CH4_conc"]
variable_attributes:
# Boudreau (1996) test cases
# Shelf/slope Rise Rise
O2_conc:initial_value: [0.250, 0.180, 0.180] # mol m-3
SO4_conc:initial_value: 28756.0e-3 # concentration mol m-3 ~ 28e-3 mol/kg * 1027 kg m-3
SO4_conc:initial_delta: 0.0
H2S_conc:initial_value: 0.0 # mol m-3
H2S_conc:initial_delta: 0.0
CH4_conc:initial_value: 0.0 # mol m-3
P_conc:initial_value: 0.001 # mol m-3

Expand Down Expand Up @@ -122,29 +124,33 @@ sediment_Corg_O2:
volume: volume_solute
variable_attributes:
R_conc:vphase: VP_Solute
R:initial_value: 0.0 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw)
R:initial_value: 1e-30 # 0.0 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw)
R:norm_value: 200e-3
R_conc:diffusivity_speciesname: O2

reservoir_SO4:
class: ReactionReservoirTotal
class: ReactionReservoirTotal
parameters:
field_data: external%SIsotope
variable_links:
R*: SO4*
volume: volume_solute
variable_attributes:
R_conc:vphase: VP_Solute
R:initial_value: 0.0 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw)
R:initial_value: 1e-30 # 0.0 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw)
R:norm_value: 20000e-3
R_conc:diffusivity_speciesname: SO4

reservoir_H2S:
class: ReactionReservoirTotal
class: ReactionReservoirTotal
parameters:
field_data: external%SIsotope
variable_links:
R*: H2S*
volume: volume_solute
variable_attributes:
R_conc:vphase: VP_Solute
R:initial_value: 0.0 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw)
R:initial_value: 1e-30 # 0.0 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw)
R:norm_value: 200e-3
R_conc:diffusivity_speciesname: H2S

Expand All @@ -155,7 +161,7 @@ sediment_Corg_O2:
volume: volume_solute
variable_attributes:
R_conc:vphase: VP_Solute
R:initial_value: 0.0 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw)
R:initial_value: 1e-30 # 0.0 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw)
R:norm_value: 200e-3
R_conc:diffusivity_speciesname: CH4

Expand All @@ -166,7 +172,7 @@ sediment_Corg_O2:
volume: volume_solute
variable_attributes:
R_conc:vphase: VP_Solute
R:initial_value: 0.0 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw)
R:initial_value: 1e-30 # 0.0 # concentration m-3 (1027 kg m-3 * 200e-6 mol/kg-sw)
R:norm_value: 2e-3
R_conc:diffusivity_speciesname: PO4

Expand All @@ -177,7 +183,7 @@ sediment_Corg_O2:
volume: volume_solid
variable_attributes:
R_conc:vphase: VP_Solid
R:initial_value: 0.0 # concentration m-3 solid phase
R:initial_value: 1e-30 # 0.0 # concentration m-3 solid phase
# mol m-3 = g cm-3 / (g mol-1) * m^3 cm-3
R:norm_value: 2083.3 # mol m-3 solid for 1% TOC = 2.5 * (1 wt%/100) / 12 * 1e6

Expand All @@ -188,7 +194,7 @@ sediment_Corg_O2:
volume: volume_solid
variable_attributes:
R_conc:vphase: VP_Solid
R:initial_value: 0.0 # concentration m-3 solid phase
R:initial_value: 1e-30 # 0.0 # concentration m-3 solid phase
# mol m-3 = g cm-3 / (g mol-1) * m^3 cm-3
R:norm_value: 2083.3 # mol m-3 solid for 1% TOC = 2.5 * (1 wt%/100) / 12 * 1e6

Expand Down
9 changes: 8 additions & 1 deletion examples/boudreau1996/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,4 +28,11 @@ Summary plots show oceanfloor solute fluxes and remineralization pathways:

![O2 and SO4 gradient summary figure](figures/PALEO_examples_sediment_x10_summary.svg)
###### Figure 1
*Oceanfloor solute fluxes and remineralization pathways vs oceanfloor [O2] and [SO4] concentration*
*Oceanfloor solute fluxes and remineralization pathways vs oceanfloor [O2] and [SO4] concentration*

# Sulphur isotope example

julia> include("PALEO_examples_sediment_Sisotopes.jl")

three sediment column example as `PALEO_examples_sediment.jl`, with sulphur isotopes enabled and low (1 mM) oceanfloor [SO4]
to illustrate Rayleigh fractionation within the sediment column as [SO4] becomes limiting.
24 changes: 24 additions & 0 deletions examples/plot_sediment.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,30 @@ function plot_solutes(
return nothing
end

function plot_solute_deltas(
output;
solutes=["SO4", "H2S"],
colT=[-1e12, 0.1, 1.0, 10.0, 100.0, 1000.0, 1e12],
colrange=1:3,
pager=PALEOmodel.DefaultPlotPager(),
)
for s in solutes
for icol in colrange
pager(plot(title="[$s] delta $icol", output, ["sediment.$(s).v_delta"], (tmodel=colT, column=icol), xlabel="$s delta (per mil)", swap_xy=true))
end
end

for s in solutes
for icol in colrange
pager(plot(title="$s flux Oceanfloor delta $icol", output, ["fluxOceanfloor.soluteflux_$s.v_delta"], (cell=icol,), ylabel="$s delta (per mil)",))
end
end

pager(:newpage)

return nothing
end

function plot_solids(
output;
solids=["Corg1"],
Expand Down

0 comments on commit 87d34fb

Please sign in to comment.