PALEOtoolkit · sjdaines · Mar 18, 2024 · Mar 18, 2024
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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