From 3f98ef5f6ad62e5ee041e52937d8893506ab00e7 Mon Sep 17 00:00:00 2001 From: Daines Date: Mon, 18 Mar 2024 19:34:19 +0000 Subject: [PATCH] Add examples for PALEOaqchem configurable chemistry --- docs/Project.toml | 6 +- docs/src/paleo_references.bib | 2 +- examples/Project.toml | 6 +- examples/configurable_chemistry/README.md | 156 ++++++++++++++++++ .../configurable_chemistry/config_ex1.yaml | 139 ++++++++++++++++ .../configurable_chemistry/config_ex2.yaml | 137 +++++++++++++++ examples/configurable_chemistry/run_ex1.jl | 72 ++++++++ examples/configurable_chemistry/run_ex2.jl | 73 ++++++++ examples/doc_order.txt | 1 + 9 files changed, 587 insertions(+), 5 deletions(-) create mode 100644 examples/configurable_chemistry/README.md create mode 100644 examples/configurable_chemistry/config_ex1.yaml create mode 100644 examples/configurable_chemistry/config_ex2.yaml create mode 100644 examples/configurable_chemistry/run_ex1.jl create mode 100644 examples/configurable_chemistry/run_ex2.jl diff --git a/docs/Project.toml b/docs/Project.toml index 1583178..0f95a0e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,6 +2,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +PALEOaqchem = "673cec3b-17d1-411f-9fcd-71c01c593120" PALEOboxes = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4" PALEOmodel = "bf7b4fbe-ccb1-42c5-83c2-e6e9378b660c" PALEOocean = "41de04b1-2efd-44ae-92ae-39d71a4fd99b" @@ -12,8 +13,9 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" [compat] Documenter = "1" DocumenterCitations = "1" -PALEOboxes = "0.20.4, 0.21" -PALEOmodel = "0.15.3" +PALEOboxes = "0.21.23" +PALEOmodel = "0.15.40" +PALEOaqchem = "0.3.9" Revise = "3.1" julia = "1.6" diff --git a/docs/src/paleo_references.bib b/docs/src/paleo_references.bib index 081397c..1de8fc9 100644 --- a/docs/src/paleo_references.bib +++ b/docs/src/paleo_references.bib @@ -35,7 +35,7 @@ @book{Zeebe2001 @article{Zhang2020, -author = {Zhang, Feifei and Lenton, Timothy M and del Rey, {\'{A}}lvaro and Romaniello, Stephen J. and Chen, Xinming and Planavsky, Noah J. and Clarkson, Matthew O. and Dahl, Tais W. and Lau, Kimberly V. and Wang, Wenqian and Li, Ziheng and Zhao, Mingyu and Isson, Terry and Algeo, Thomas J. and Anbar, Ariel D.}, +author = {Zhang, Feifei and Lenton, Timothy M and del Rey, Alvaro and Romaniello, Stephen J. and Chen, Xinming and Planavsky, Noah J. and Clarkson, Matthew O. and Dahl, Tais W. and Lau, Kimberly V. and Wang, Wenqian and Li, Ziheng and Zhao, Mingyu and Isson, Terry and Algeo, Thomas J. and Anbar, Ariel D.}, doi = {10.1016/j.gca.2020.05.011}, journal = {Geochimica et Cosmochimica Acta}, month = {oct}, diff --git a/examples/Project.toml b/examples/Project.toml index 2022211..66b355f 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -5,6 +5,7 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" NBInclude = "0db19996-df87-5ea3-a455-e3a50d440464" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +PALEOaqchem = "673cec3b-17d1-411f-9fcd-71c01c593120" PALEOboxes = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4" PALEOcopse = "4a6ed817-0e58-48c6-8452-9e9afc8cb508" PALEOmodel = "bf7b4fbe-ccb1-42c5-83c2-e6e9378b660c" @@ -15,5 +16,6 @@ Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" [compat] -PALEOboxes = "0.20.4, 0.21" -PALEOmodel = "0.15.3" +PALEOboxes = "0.21.23" +PALEOmodel = "0.15.40" +PALEOaqchem = "0.3.9" diff --git a/examples/configurable_chemistry/README.md b/examples/configurable_chemistry/README.md new file mode 100644 index 0000000..a53eaed --- /dev/null +++ b/examples/configurable_chemistry/README.md @@ -0,0 +1,156 @@ +# Configurable chemistry using PALEOaqchem + +These examples illustrate how to use the yaml configuration file to configure generic chemistry Reactions from PALEOaqchem +to implement the same minimal model of the marine carbonate system shown in the [Ocean chemistry](@ref) examples. + +See `PALEOaqchem` [documentation](https://paleotoolkit.github.io/PALEOaqchem.jl/) for more detail on the generic aqueous chemistry. + +See `PALEOmodel` [documentation](https://paleotoolkit.github.io/PALEOmodel.jl/) for more information on analysing model output. + +## Example 1 Configuring a minimal Alk-pH model + +This example configures generic chemistry reactions from PALEOaqchem to define a minimal Alk-pH model for aqueous carbonate chemistry. Carbonate system equations are from [Zeebe2001](@cite). + +The `ocean` domain contains a [`ReactionNoTransport`](https://paleotoolkit.github.io/PALEOocean.jl/dev/PALEOocean_Reactions/#PALEOocean.Ocean.OceanNoTransport.ReactionOceanNoTransport) from the [`PALEOocean`](https://github.com/PALEOtoolkit/PALEOocean.jl) Julia package. This defines standard ocean variables including cell `volume`, and is configured here to provide one ocean cell. + +There are two state variables for `DIC` and `TAlk`, implemented using a [`ReactionReservoir`](https://paleotoolkit.github.io/PALEOboxes.jl/stable/ReactionCatalog/#PALEOboxes.Reservoirs.ReactionReservoir) from the [`PALEOboxes`](https://github.com/PALEOtoolkit/PALEOboxes.jl) Julia package. This Reaction also provides concentrations in `mol m-3`. Constant Boron concentration `B_conc` is defined using a [`ReactionReservoirConst`](https://paleotoolkit.github.io/PALEOboxes.jl/stable/ReactionCatalog/#PALEOboxes.Reservoirs.ReactionReservoirConst). + +The carbonate chemistry system includes the speciation of water, boron and DIC. `B`, `DIC`, and `TAlk` are the three total variables with corresponding primary species concentrations `BOH3_conc`, `CO2_aq_conc`, and `H_conc` (defined by `pH`). Primary species are defined using the PALEOaqchem `ReactionConstraintReservoir`, which defines a state variable (attribute `:vfunction = PB.VF_State`) eg `BOH3_conc`, an algebraic constraint for the corresponding total eg `B_constraint_conc` (with attribute `:vfunction = PB.VF_Constraint`), and a variable eg `B_calc` to accumulate primary and secondary species contributions. Secondary species eg `BOH4` are defined using the PALEOaqchem `ReactionAqEqb`, and their contribution added to the calculated total eg `B_calc`. The PALEO DAE solver then solves for the primary species concentration `BOH3_conc` that makes `B_calc` equal to the required value `B`. Note that the `ReactionConstraintReservoir` for `TAlk` defines `pH` as the state variable (from which `H_conc` is calculated), and that many species contribute to both `TAlk_calc` and their own total. + +The combined system is then a Differential Algebraic Equation (DAE) system with two ODE state variables and corresponding time derivatives for `DIC` and `TAlk`, three algebraic constraints for `B`, `DIC`, and `TAlk`, and three additional state variables for the primary species `BOH3_conc`, `CO2_aq_conc`, and `pH`, and is integrated forward in time by a DAE solver [`PALEOmodel.ODE.integrateDAE`](https://paleotoolkit.github.io/PALEOmodel.jl/stable/PALEOmodelSolvers/#PALEOmodel.ODE.integrateDAE). + + +### yaml configuration file +The model configuration (file `examples/configurable_chemistry/config_ex1.yaml`) contains two Reservoirs `DIC`, `TAlk`. + +A `ReactionFluxPerturb` from the PALEOboxes.jl reaction catalog is used to add a constant TAlk flux. + +```@eval +import Markdown +Markdown.parse( + """```julia + $(read(joinpath(ENV["PALEO_EXAMPLES"], "configurable_chemistry/config_ex1.yaml"), String)) + ```""" +) +``` + +### Run script +The script to run the model (file `examples/configurable_chemistry/run_ex1.jl`) contains: +```@eval +import Markdown +Markdown.parse( + """```julia + $(read(joinpath(ENV["PALEO_EXAMPLES"], "configurable_chemistry/run_ex1.jl"), String)) + ```""" +) +``` +and produces output showing the change, if `TAlk_conc` increase, how the carbonic acid and pH change in ocean: +```@example ex1 +include(joinpath(ENV["PALEO_EXAMPLES"], "configurable_chemistry/run_ex1.jl")) # hide +plot(paleorun.output, ["ocean.TAlk_conc", "ocean.DIC_conc"], (cell=1,); ylabel="TAlk, DIC conc (mol m-3)") # hide +savefig("ex1_plot1.svg"); nothing # hide +plot(paleorun.output, "ocean.pH", (cell=1,)) # hide +savefig("ex1_plot2.svg"); nothing # hide +plot(paleorun.output, ["ocean.DIC_conc", "ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc"], (cell=1,); ylabel="DIC species (mol m-3)") # hide +savefig("ex1_plot3.svg"); nothing # hide +display( # hide + plot( # hide + paleorun.output, # hide + ["ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc", "ocean.BOH4_conc", "ocean.BOH3_conc", "ocean.H_conc", "ocean.OH_conc",], # hide + (cell=1, tmodel=(1.0, 1e12)); # omit first point (pH 8 starting condition) # hide + coords=["tmodel"=>("ocean.pH",),], # plot against pH instead of tmodel # hide + ylabel="H2O, B, DIC species (mol m-3)", ylim=(0.5e-3, 0.5e1), yscale=:log10, # hide + legend_background_color=nothing, # hide + legend=:bottom, # hide + ) # hide +) # hide +savefig("ex1_plot4.svg"); nothing # hide +``` + +![](ex1_plot1.svg) +![](ex1_plot2.svg) +![](ex1_plot3.svg) +![](ex1_plot4.svg) + + +### Displaying model structure and variables + +All metadata for model variables can be displayed with `PB.show_variables`: +```@example ex1 +show(PB.show_variables(model), allcols=true, allrows=true) # display in REPL +# vscodedisplay(PB.show_variables(model)) # more convenient when using VS code +``` + +## Example 2 Minimal Alk-pH model with implicit variables + +This example contains the same chemistry as [Example 1 Configuring a minimal Alk-pH model](@ref), using the Direct Substitution Approach (DSA) to eliminate the two algebraic equations corresponding to the mass action equations for `DIC` and `TAlk`. + +The state variables for `DIC_conc` and `TAlk_conc` are now implicit total variables defined by PALEOaqchem `ReactionImplicitReservoir`s. These define a PALEO Total variable (attribute `:vfunction = PB.VF_Total`) eg `DIC_conc` which is a function of a primary species concentration implemented as a PALEO StateTotal variable (attribute `:vfunction = PB.VF_StateTotal`) eg `CO_2_aq_conc`, and a time derivative eg `DIC_conc_sms` implemented as a PALEO Deriv variable (attribute `:vfunction = PB.VF_Deriv`). + +The combined system is then a Differential Algebraic Equation (DAE) system with three implicit total variables for `DIC_conc`, `B_conc`, and `TAlk_conc`, with corresponding time derivatives `DIC_conc_sms`, `B_conc_sms`, and `TAlk_conc_sms` and state variables for the primary species `CO2_aq_conc`, `BOH3_conc` and `pH`, and is integrated forward in time by a DAE solver [`PALEOmodel.ODE.integrateDAE`](https://paleotoolkit.github.io/PALEOmodel.jl/stable/PALEOmodelSolvers/#PALEOmodel.ODE.integrateDAE). + +Notes: +- The combination of constant `BOH3_conc` and a `ReactionConstraintReservoir` defining primary species `BOH3_conc` has been replaced with an implicit state variable for `B_conc` and primary species `BOH3_conc`: this is required to use the IDA DAE solver, which cannot handle initialisation of that combination of implicit and constraint variables (see `PALEOmodel` documentation). +- Initial values are now specified for the primary species (eg `CO2_conc`, `BOH3_conc`, `pH`), not the corresponding totals (eg `DIC_conc`). This is less convenient here (where we are want to set `DIC_conc` and `B_conc` that remain constant), but may not be an issue in a more complete model that includes additional processes that define time evolution that is insensitive to initial conditions after a spinup. +- This configuration provides concentrations, instead of mol (per cell) to the solver. This is not required here, but is good practice as it helps the numerical solver by providing variables with better scaling. + + +### yaml configuration file +The model configuration (file `examples/configurable_chemistry/config_ex2.yaml`) contains: + +```@eval +import Markdown +Markdown.parse( + """```julia + $(read(joinpath(ENV["PALEO_EXAMPLES"], "configurable_chemistry/config_ex2.yaml"), String)) + ```""" +) +``` + +### Run script +The script to run the model (file `examples/configurable_chemistry/run_ex2.jl`) contains: +```@eval +import Markdown +Markdown.parse( + """```julia + $(read(joinpath(ENV["PALEO_EXAMPLES"], "configurable_chemistry/run_ex2.jl"), String)) + ```""" +) +``` +and produces output showing the change, if `TAlk_conc` increase, how the carbonic acid and pH change in ocean: +```@example ex2 +include(joinpath(ENV["PALEO_EXAMPLES"], "configurable_chemistry/run_ex2.jl")) # hide +plot(paleorun.output, ["ocean.TAlk_conc", "ocean.DIC_conc"], (cell=1,); ylabel="TAlk, DIC conc (mol m-3)") # hide +savefig("ex2_plot1.svg"); nothing # hide +plot(paleorun.output, "ocean.pH", (cell=1,)) # hide +savefig("ex2_plot2.svg"); nothing # hide +plot(paleorun.output, ["ocean.DIC_conc", "ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc"], (cell=1,); ylabel="DIC species (mol m-3)") # hide +savefig("ex2_plot3.svg"); nothing # hide +display( # hide + plot( # hide + paleorun.output, # hide + ["ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc", "ocean.BOH4_conc", "ocean.BOH3_conc", "ocean.H_conc", "ocean.OH_conc",], # hide + (cell=1, tmodel=(1.0, 1e12)); # omit first point (pH 8 starting condition) # hide + coords=["tmodel"=>("ocean.pH",),], # plot against pH instead of tmodel # hide + ylabel="H2O, B, DIC species (mol m-3)", ylim=(0.5e-3, 0.5e1), yscale=:log10, # hide + legend_background_color=nothing, # hide + legend=:bottom, # hide + ) # hide +) # hide +savefig("ex2_plot4.svg"); nothing # hide +``` + +![](ex2_plot1.svg) +![](ex2_plot2.svg) +![](ex2_plot3.svg) +![](ex2_plot4.svg) + + +### Displaying model structure and variables + +All metadata for model variables can be displayed with `PB.show_variables`: +```@example ex2 +show(PB.show_variables(model), allcols=true, allrows=true) # display in REPL +# vscodedisplay(PB.show_variables(model)) # more convenient when using VS code +``` + diff --git a/examples/configurable_chemistry/config_ex1.yaml b/examples/configurable_chemistry/config_ex1.yaml new file mode 100644 index 0000000..4c285ad --- /dev/null +++ b/examples/configurable_chemistry/config_ex1.yaml @@ -0,0 +1,139 @@ +Minimal_Alk_pH: + + domains: + + global: + + reactions: + add_Alk: + # from PALEOboxes reaction catalog: constant input flux of Alk + class: ReactionFluxPerturb + parameters: + # linear interpolation for input flux - set to constant + perturb_times: [-1e30, 1e30] + # ocean volume = 3.6e14 * 3688.0 = 1.33e18 m^3 + # DIC in ocean = 2 * 1.33e18 = 2.66e18 mol + # set Alk flux so Alk = DIC after 100 yr-1 = 2.66e18/100 = 2.66e16 mol yr-1 + perturb_totals: [2.66e16, 2.66e16] + variable_links: + F: ocean.TAlk_sms # add flux to TAlk reservoir NB: works for single-box ocean only !! + # -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + + ocean: + + reactions: + + oceantrasport: + + class: ReactionOceanNoTransport + + parameters: + + area: [3.6e14] + depth: [3688.0] + + + + reservoir_TA: + + class: ReactionReservoir + + variable_links: + R*: TAlk* + variable_attributes: + R:initial_value: 0.0 + R:norm_value: 2.4e0 + + reservoir_DIC: + class: ReactionReservoir + + variable_links: + R*: DIC* + variable_attributes: + R:initial_value: 2.0e0 + R:norm_value: 2.0e0 + + constant_B: + class: ReactionReservoirConst + variable_links: + R*: B + variable_attributes: + R_conc:initial_value: 0.427 # mol m-3 contemporary value + + ########################################################################## + # TAlk, primary species H_conc (as pH) + ######################################################################## + + # [H+] primary species (as pH) for TAlk total + H_primary_species: {class: ReactionConstraintReservoir, + variable_links: {Primary_pconc: pH, Primary_conc: H_conc, R*: TAlk*}, + parameters: {primary_total_stoich: -1.0, primary_variable: p_concentration, constraint_variable: concentration}, + variable_attributes: {Primary_pconc%initial_value: 7.0, Primary_pconc%norm_value: 1.0, R_constraint_conc%norm_value: 1e-3}} + + # Define OH_conc and add to TAlk + # OH- + H+ <--> H2O, [OH] = K_w / [H] + OH_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator + K_eqb: 6.0e-14, # K_w 6.0e-14 (mol^2 kg-2) equilibrium constant of water at S=35, T=25°C, + K_density_power: 2.0, # K_eqb*rho_ref^2 to convert to (mol^2 m-6) + Reactants: ["OH_conc", "H_conc"], Products: [], + N_components: ["TAlk_calc"], + } + } + + ########################################################################## + # Boron speciation + ########################################################################## + + # BOH3 primary species for B total + BOH3_primary_species: {class: ReactionConstraintReservoir, + variable_links: {Primary*: BOH3*, R*: B*}, + parameters: {primary_total_stoich: 1.0, primary_variable: concentration, constraint_variable: concentration}, + variable_attributes: {Primary_conc%initial_value: 1.0, Primary_conc%norm_value: 1e-3, R_constraint_conc%norm_value: 1e-3}} + + + # Define BOH4_conc and add to B, TAlk + # B(OH)_4^- + H+ <--> B(OH)3 + H2O, [B(OH)_4^-] = K_B [B(OH)3] / [H+] + BOH4_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator + K_eqb: 2.5e-9, # (mol kg-1) K_B equilibrium constant of B(OH)4- + K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3) + Reactants: ["BOH4_conc", "H_conc"], Products: ["BOH3_conc"], + N_components: ["B_calc", "TAlk_calc"], + } + } + + ########################################################################## + # Carbon speciation + ########################################################################## + + # CO2_aq primary species for DIC total + CO2_aq_primary_species: {class: ReactionConstraintReservoir, + variable_links: {Primary*: CO2_aq*, R*: DIC*}, + parameters: {primary_total_stoich: 1.0, primary_variable: concentration, constraint_variable: concentration}, + variable_attributes: {Primary_conc%initial_value: 1.0, Primary_conc%norm_value: 1e-3, R_constraint_conc%norm_value: 1.0}} + + + # Define HCO3_conc and add to DIC, TAlk + # HCO3- + H+ <--> CO2 + H2O, [HCO3-] = K_1 [CO2] / [H+] + HCO3_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator + K_eqb: 1.4e-6, # (mol kg-1) K_1 equilibrium constant of CO2_aq and HCO3- + K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3) + Reactants: ["HCO3_conc", "H_conc"], Products: ["CO2_aq_conc"], + N_components: ["DIC_calc", "TAlk_calc"], + } + } + + # Define CO3_conc and add to DIC, TAlk + # CO3-- + H+ <--> HCO3- + H2O, [CO3--] = K_2 [HCO3-] / [H+] + CO3_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator + K_eqb: 1.2e-9, # (mol kg-1) K_2 equilibrium constant of HCO3- and CO32- + K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3) + Reactants: ["CO3_conc", "H_conc"], Products: ["HCO3_conc"], + N_components: ["DIC_calc", "2*TAlk_calc"], + } + } + + oceanfloor: + # unused here, set up by ReactionOceanNoTransport + + oceansurface: + # unused here, set up by ReactionOceanNoTransport \ No newline at end of file diff --git a/examples/configurable_chemistry/config_ex2.yaml b/examples/configurable_chemistry/config_ex2.yaml new file mode 100644 index 0000000..a0faf9b --- /dev/null +++ b/examples/configurable_chemistry/config_ex2.yaml @@ -0,0 +1,137 @@ +Minimal_Alk_pH: + + domains: + + global: + + reactions: + add_Alk: + # from PALEOboxes reaction catalog: constant input flux of Alk + class: ReactionFluxPerturb + parameters: + # linear interpolation for input flux - set to constant + perturb_times: [-1e30, 1e30] + # ocean volume = 3.6e14 * 3688.0 = 1.33e18 m^3 + # DIC in ocean = 2 * 1.33e18 = 2.66e18 mol + # set Alk flux so Alk = DIC after 100 yr-1 = 2.66e18/100 = 2.66e16 mol yr-1 + perturb_totals: [2.66e16, 2.66e16] + variable_links: + F: ocean.TAlk_sms # add flux to TAlk reservoir NB: works for single-box ocean only !! + # -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + + ocean: + + reactions: + + oceantransport: + + class: ReactionOceanNoTransport + + parameters: + + area: [3.6e14] + depth: [3688.0] + + # TA as total and H (as pH) as primary species + TA_total_H_primary: + class: ReactionImplicitReservoir + variable_links: + Primary_pconc: pH + Primary_conc: H_conc + R*: TAlk* + parameters: + primary_total_stoich: -1.0 + primary_variable: p_concentration # provide solver with -log10(H_conc) + total_variable: concentration # provide solver with TAlk_conc + variable_attributes: + Primary_pconc%initial_value: 4.28 # value tuned to give TAlk ~ 0.0 initially + Primary_pconc%norm_value: 1.0 + R_conc%norm_value: 1.0 + + # DIC as total and CO2_aq as primary species + DIC_total_CO2_aq_primary: + class: ReactionImplicitReservoir + variable_links: + R*: DIC* + Primary*: CO2_aq* + parameters: + primary_total_stoich: 1.0 # add CO2 to DIC, no contribution to TAlk + primary_variable: concentration # provide solver with CO2_aq_conc (mol m-3) + total_variable: concentration # provide solver with DIC_conc (mol m-3) + variable_attributes: + Primary_conc%initial_value: 1.94 # This value has been tuned to produce initial DIC_conc = 2.0 + Primary_conc%norm_value: 1.0 + R_conc%norm_value: 1.0 + + # B as total and BOH3 as primary species + B_total_BOH3_primary: + class: ReactionImplicitReservoir + variable_links: + R*: B* + Primary*: BOH3* + parameters: + primary_total_stoich: 1.0 # add BOH3 to B, no contribution to TAlk + primary_variable: concentration # provide solver with BOH3_conc (mol m-3) + total_variable: concentration # provide solver with B_conc (mol m-3) + variable_attributes: + Primary_conc%initial_value: 0.427 # This produces initial B_conc ~ 0.427 (as initial BOH4_conc is very small) + Primary_conc%norm_value: 1.0 + R_conc%norm_value: 1.0 + + ########################################################################## + # H2O speciation + ######################################################################## + + # Define OH_conc and add to TAlk + # OH- + H+ <--> H2O, [OH] = K_w / [H] + OH_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator + K_eqb: 6.0e-14, # K_w 6.0e-14 (mol^2 kg-2) equilibrium constant of water at S=35, T=25°C, + K_density_power: 2.0, # K_eqb*rho_ref^2 to convert to (mol^2 m-6) + Reactants: ["OH_conc", "H_conc"], Products: [], + N_components: ["TAlk_calc"], + } + } + + ########################################################################## + # Boron speciation + ########################################################################## + + # Define BOH4_conc and add to B, TAlk + # B(OH)_4^- + H+ <--> B(OH)3 + H2O, [B(OH)_4^-] = K_B [B(OH)3] / [H+] + BOH4_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator + K_eqb: 2.5e-9, # (mol kg-1) K_B equilibrium constant of B(OH)4- + K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3) + Reactants: ["BOH4_conc", "H_conc"], Products: ["BOH3_conc"], + N_components: ["B_calc", "TAlk_calc"], + } + } + + ########################################################################## + # Carbon speciation + ########################################################################## + + # Define HCO3_conc and add to DIC, TAlk + # HCO3- + H+ <--> CO2 + H2O, [HCO3-] = K_1 [CO2] / [H+] + HCO3_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator + K_eqb: 1.4e-6, # (mol kg-1) K_1 equilibrium constant of CO2_aq and HCO3- + K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3) + Reactants: ["HCO3_conc", "H_conc"], Products: ["CO2_aq_conc"], + N_components: ["DIC_calc", "TAlk_calc"], + } + } + + # Define CO3_conc and add to DIC, TAlk + # CO3-- + H+ <--> HCO3- + H2O, [CO3--] = K_2 [HCO3-] / [H+] + CO3_conc: { class: ReactionAqEqb, parameters: { K_power: 1.0, # we need K_eqb in numerator + K_eqb: 1.2e-9, # (mol kg-1) K_2 equilibrium constant of HCO3- and CO32- + K_density_power: 1.0, # K_eqb*rho_ref to convert to (mol m-3) + Reactants: ["CO3_conc", "H_conc"], Products: ["HCO3_conc"], + N_components: ["DIC_calc", "2*TAlk_calc"], + } + } + + oceanfloor: + # unused here, set up by ReactionOceanNoTransport + + oceansurface: + # unused here, set up by ReactionOceanNoTransport \ No newline at end of file diff --git a/examples/configurable_chemistry/run_ex1.jl b/examples/configurable_chemistry/run_ex1.jl new file mode 100644 index 0000000..03c3b2b --- /dev/null +++ b/examples/configurable_chemistry/run_ex1.jl @@ -0,0 +1,72 @@ +import PALEOboxes as PB +import PALEOmodel + +import PALEOocean +import PALEOaqchem +using Plots + + + + +##################################################### +# Create model + +model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex1.yaml"), "Minimal_Alk_pH") + +######################################################### +# Initialize +########################################################## + +initial_state, modeldata = PALEOmodel.initialize!(model) + +##################################################################### +# Optional: call ODE function to check derivative +####################################################################### +initial_deriv = similar(initial_state) +PALEOmodel.ODE.ModelODE(modeldata)(initial_deriv, initial_state , nothing, 0.0) +println("initial_state: ", initial_state) +println("initial_deriv: ", initial_deriv) + +################################################################# +# Integrate vs time +################################################################## + +# create a Run object to hold model and output +paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + +println("integrate, DAE") +# first run is slow as it includes JIT time +@time PALEOmodel.ODE.integrateDAE( + paleorun, initial_state, modeldata, (0.0, 200.0), + solvekwargs=( + reltol=1e-6, + # saveat=0.1, # save output every 0.1 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/ + ) +); + +######################################## +# Plot output +######################################## + +display(plot(paleorun.output, ["ocean.TAlk_conc", "ocean.DIC_conc"], (cell=1,); + ylabel="TAlk, DIC conc (mol m-3)")) +display(plot(paleorun.output, ["ocean.TAlk","ocean.TAlk_sms"], (cell=1,))) +display(plot(paleorun.output, "ocean.pH", (cell=1,))) +display(plot(paleorun.output, ["ocean.H_conc", "ocean.OH_conc"], (cell=1,); + ylabel="H2O species (mol m-3)", ylim=(0, Inf))) +display(plot(paleorun.output, ["ocean.B_conc", "ocean.BOH4_conc", "ocean.BOH3_conc"], (cell=1,); + ylabel="B species (mol m-3)", ylim=(0, Inf))) +display(plot(paleorun.output, ["ocean.DIC_conc", "ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc"], (cell=1,); + ylabel="DIC species (mol m-3)", ylim=(0, Inf))) +display( + plot( + paleorun.output, + ["ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc", "ocean.BOH4_conc", "ocean.BOH3_conc", "ocean.H_conc", "ocean.OH_conc",], + (cell=1, tmodel=(1.0, 1e12)); # omit first point (pH 8 starting condition) + coords=["tmodel"=>("ocean.pH",),], # plot against pH instead of tmodel + ylabel="H2O, B, DIC species (mol m-3)", ylim=(0.5e-3, 0.5e1), yscale=:log10, + legend_background_color=nothing, + legend=:bottom, + ) +) + diff --git a/examples/configurable_chemistry/run_ex2.jl b/examples/configurable_chemistry/run_ex2.jl new file mode 100644 index 0000000..56dc89f --- /dev/null +++ b/examples/configurable_chemistry/run_ex2.jl @@ -0,0 +1,73 @@ +import PALEOboxes as PB +import PALEOmodel + +import PALEOocean +import PALEOaqchem +using Plots + + + + +##################################################### +# Create model + +model = PB.create_model_from_config(joinpath(@__DIR__, "config_ex2.yaml"), "Minimal_Alk_pH") + +######################################################### +# Initialize +########################################################## + +initial_state, modeldata = PALEOmodel.initialize!(model) + +##################################################################### +# Optional: call ODE function to check derivative +####################################################################### +initial_deriv = similar(initial_state) +PALEOmodel.ODE.ModelODE(modeldata)(initial_deriv, initial_state , nothing, 0.0) +println("initial_state: ", initial_state) +println("initial_deriv: ", initial_deriv) + +################################################################# +# Integrate vs time +################################################################## + +# create a Run object to hold model and output +paleorun = PALEOmodel.Run(model=model, output = PALEOmodel.OutputWriters.OutputMemory()) + +println("integrate, DAE") +# first run is slow as it includes JIT time +# NB: Jacobian is required to use implicit variables, so use integrateDAEForwardDiff instead of integrateDAE +@time PALEOmodel.ODE.integrateDAEForwardDiff( + paleorun, initial_state, modeldata, (0.0, 200.0), + solvekwargs=( + reltol=1e-6, + # saveat=0.1, # save output every 0.1 yr see https://diffeq.sciml.ai/dev/basics/common_solver_opts/ + ) +); + +######################################## +# Plot output +######################################## + +display(plot(paleorun.output, ["ocean.TAlk_conc", "ocean.DIC_conc"], (cell=1,); + ylabel="TAlk, DIC conc (mol m-3)")) +display(plot(paleorun.output, ["ocean.TAlk_conc","ocean.TAlk_conc_sms"], (cell=1,))) +display(plot(paleorun.output, "ocean.pH", (cell=1,))) +display(plot(paleorun.output, ["ocean.H_conc", "ocean.OH_conc"], (cell=1,); + ylabel="H2O species (mol m-3)", ylim=(0, Inf))) +display(plot(paleorun.output, ["ocean.B_conc", "ocean.BOH4_conc", "ocean.BOH3_conc"], (cell=1,); + ylabel="B species (mol m-3)", ylim=(0, Inf))) +display(plot(paleorun.output, ["ocean.DIC_conc", "ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc"], (cell=1,); + ylabel="DIC species (mol m-3)", ylim=(0, Inf))) +display( + plot( + paleorun.output, + ["ocean.HCO3_conc", "ocean.CO3_conc", "ocean.CO2_aq_conc", "ocean.BOH4_conc", "ocean.BOH3_conc", "ocean.H_conc", "ocean.OH_conc",], + (cell=1, tmodel=(1.0, 1e12)); # omit first point (pH 8 starting condition) + coords=["tmodel"=>("ocean.pH",),], # plot against pH instead of tmodel + ylabel="H2O, B, DIC species (mol m-3)", ylim=(0.5e-3, 0.5e1), yscale=:log10, + legend_background_color=nothing, + legend=:bottom, + ) +) + diff --git a/examples/doc_order.txt b/examples/doc_order.txt index cdf1894..972a6bc 100644 --- a/examples/doc_order.txt +++ b/examples/doc_order.txt @@ -4,4 +4,5 @@ error_examples CPU CPU_modular ocean_chemistry +configurable_chemistry