Skip to content

Commit

Permalink
Add some more doc for Xiyuan's ocean chemistry examples
Browse files Browse the repository at this point in the history 
- more words to explain how PALEO works
- generalize air-sea to gas X (not just CO2!)
  • Loading branch information
@sjdaines
sjdaines committed May 16, 2023
1 parent d1bc22f commit 9268449
Show file tree
Hide file tree
Showing 8 changed files with 138 additions and 81 deletions.
11 changes: 11 additions & 0 deletions docs/src/paleo_references.bib
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,17 @@ @article{Clarkson2018
pages = {2918--2923},
}


@book{Zeebe2001,
title = {{CO2} in {Seawater}: {Equilibrium}, {Kinetics}, {Isotopes}},
isbn = {978-0-444-50946-8},
url = {https://www.elsevier.com/books/co2-in-seawater-equilibrium-kinetics-isotopes-65/zeebe/978-0-444-50946-8},
publisher = {Elsevier},
author = {Zeebe, Richard E and Wolf-Gladrow, Dieter a.},
year = {2001},
}


@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.},
doi = {10.1016/j.gca.2020.05.011},
Expand Down
1 change: 1 addition & 0 deletions examples/doc_order.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,5 +2,6 @@ reservoirs
solvers
error_examples
CPU
CPU_modular
ocean_chemistry

45 changes: 41 additions & 4 deletions examples/ocean_chemistry/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,10 +8,16 @@ See `PALEOmodel` [documentation](https://paleotoolkit.github.io/PALEOmodel.jl/)

## Example 1 Minimal Alk-pH model

Generalizing the Reaction to establish a minimal Alk-pH model. Math equations are from Richard E. Zeebe's book(2001).
This example implements a new `Reaction_Alk_pH` to establish 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`.

The carbonate chemistry system is solved as an algebraic constraint, calculating `TAlk_calcu` as a function of `pH` and then using a PALEO solver to solve for the `pH` value that makes `TAlk_calcu` equal to the required value. In PALEO this is implemented by defining `pH` as a `VarState` (attribute `:vfunction = PB.VF_State`) and the algebraic constraint `TAlk_error` as a `VarConstraint` (with attribute `:vfunction = PB.VF_Constraint`). The combined system of `TAlk` and `DIC` reservoirs and `pH` is then a Differential Algebraic Equation (DAE) with two state variables and one algebraic constraint, and is integrated forward in time by a DAE solver []`PALEOmodel.ODE.integrateDAE`](https://paleotoolkit.github.io/PALEOmodel.jl/stable/PALEOmodelSolvers/#PALEOmodel.ODE.integrateDAE).

### Additional code files
The Reaction code (file `examples/ocean_chemistry/reactions_Alk_pH.jl`) now produces calculation of different carbonic acid `HCO3_conc`, `CO3_conc`, `CO2_aq_conc`, and coeval concentrations change of `BOH4_conc`, `H_conc`, `OH_conc`, `DIC_conc` and `TAlk_conc`:
The Reaction code (`Reaction_Alk_pH` in file `examples/ocean_chemistry/reactions_Alk_pH.jl`) now produces calculation of different carbon species `HCO3_conc`, `CO3_conc`, `CO2_aq_conc`, boron species `BOH4_conc`, water species `H_conc` and `OH_conc` given `DIC_conc`, `TAlk_conc` and `pH`. Difference from required alkalinity `TAlk_error` is then calculated and labelled as an algebraic constraint. Note that there is loop over ocean cells to support arbitrary ocean models:
```@eval
str = read("../../../../examples/ocean_chemistry/reactions_Alk_pH.jl", String)
str = """```julia
Expand All @@ -21,6 +27,14 @@ import Markdown
Markdown.parse(str)
```

Documentation (generated by the Julia docstring) reads:
```@meta
CurrentModule = Min_Alk_pH
```
```@docs
Reaction_Alk_pH
```

### yaml configuration file
The model configuration (file `examples/ocean_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.
Expand Down Expand Up @@ -52,11 +66,24 @@ plot(paleorun.output, "ocean.pH",
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
Expand All @@ -69,11 +96,13 @@ show(PB.show_variables(model), allcols=true, allrows=true) # display in REPL

## Example 2 Air-sea exchange

Adds air-sea exchange of CO2 to Example 1
This example adds air-sea exchange of CO2 to [Example 1 Minimal Alk-pH model](@ref).

This configuration adds an atmosphere Domain `atm` with a state variable for atmospheric CO2. Following the standard PALEO convention for [`coupling spatial Domains`](https://paleotoolkit.github.io/PALEOboxes.jl/stable/DesignOverview/#Coupling-Spatial-Domains), air-sea exchange is implemented by a combination of the new reaction `Reaction_Min_AirSeaExchange` in the `oceansurface` Domain, a [`ReactionFluxTarget`](https://paleotoolkit.github.io/PALEOboxes.jl/stable/ReactionCatalog/#PALEOboxes.Fluxes.ReactionFluxTarget) in a `fluxAtmtoOceansurface` Domain to store the calculated flux, and a pair of [`ReactionFluxTransfer`s](https://paleotoolkit.github.io/PALEOboxes.jl/stable/ReactionCatalog/#PALEOboxes.Fluxes.ReactionFluxTransfer) to apply the calculated fluxes to the atmosphere CO2 and ocean DIC reservoirs. NB: the `ocean.oceansurface` subdomain represents the subset of ocean cells adjacent to the surface, and contains the same number of cells as the `oceansurface` and `fluxAtmtoOceansurface` Domains, with a 1-1 correspondence.

### Additional code files

In order to evaluate the CO2 flux change between air and sea, we add a file (file `examples/ocean_chemistry/reactions_AirSeaExchange.jl`) to achieve air-sea CO2 exchange following Henry's Law.
In order to evaluate the CO2 flux change between air and sea, we add a file (file `examples/ocean_chemistry/reactions_AirSeaExchange.jl`) that implements a `Reaction_Min_AirSeaExchange` to calculate air-sea CO2 exchange following Henry's Law. NB: the reaction is implemented for a generic gas `X` that is then linked to the CO2 and DIC variables using the `variable_links:` section in the .yaml configuration file.
```@eval
str = read("../../../../examples/ocean_chemistry/reactions_AirSeaExchange.jl", String)
str = """```julia
Expand All @@ -83,6 +112,14 @@ import Markdown
Markdown.parse(str)
```

Documentation (generated by the Julia docstring) reads:
```@meta
CurrentModule = Min_AirSeaExchange
```
```@docs
Reaction_Min_AirSeaExchange
```

### yaml configuration file
The model configuration (file `examples/ocean_chemistry/config_ex2.yaml`) contains three Reservoirs `DIC`, `TAlk` and `CO2`. Following `reservoirs` [Example 4 Transfer between Domains](@ref), we use `ReactionFluxTarget` and `ReactionFluxTransfer` to transfer `CO2_airsea_exchange` between `DIC` reservoir and `CO2` reservoir:
```@eval
Expand Down
13 changes: 7 additions & 6 deletions examples/ocean_chemistry/config_ex1.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,22 +55,23 @@ Minimal_Alk_pH:
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

solve_Alk_pH:

class: Reaction_Alk_pH

parameters:

# kappa: 1.0
# density: 1027.0
# TAlk_conc_change: 10.0e-7
parameters:

K_1: 1.4e-6
K_2: 1.2e-9
K_w: 6.0e-14
K_B: 2.5e-9
B_total: 4.2e-4

variable_links:
density: rho_ref # OceanNoTransport provides density as rho_ref
Expand Down
23 changes: 12 additions & 11 deletions examples/ocean_chemistry/config_ex2.yaml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,22 +72,23 @@ Minimal_Alk_pH_AirSea:
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

solve_Alk_pH:

class: Reaction_Alk_pH

parameters:

# kappa: 1.0
# density: 1027.0
# TAlk_conc_change: 10.0e-7
parameters:

K_1: 1.4e-6
K_2: 1.2e-9
K_w: 6.0e-14
K_B: 2.5e-9
B_total: 4.2e-4
K_B: 2.5e-9

variable_links:
density: rho_ref # OceanNoTransport provides density as rho_ref
Expand All @@ -108,16 +109,16 @@ Minimal_Alk_pH_AirSea:
parameters:

K_0: 3.4e1 # (mol m-3 atm-1) = 3.4e-2 (mol L-1 atm-1) * 1e3 (L m-3)
vpiston: 1138.8
vpiston: 1138.8 # m yr-1

variable_links:

pCO2atm: atm.pCO2atm
CO2_aq_conc: ocean.oceansurface.CO2_aq_conc
pXatm: atm.pCO2atm
X_aq_conc: ocean.oceansurface.CO2_aq_conc

area: Asurf # ocean surface area

CO2_airsea_exchange: fluxAtmtoOceansurface.flux_CO2
X_airsea_exchange: fluxAtmtoOceansurface.flux_CO2


# apply air-sea flux to ocean surface
Expand Down
35 changes: 21 additions & 14 deletions examples/ocean_chemistry/reactions_AirSeaExchange.jl
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,24 @@
module Min_AirSeaExchange

import PALEOboxes as PB
using PALEOboxes.DocStrings # for $(PARS) and $(METHODS_DO)

"""
Min_AirSeaExchange
Reaction_Min_AirSeaExchange
Minimal example, just make a easy way to illustrate Air-Sea exchange.
Current this file only consider CO2 exchange between Air and Sea.
This implements exchange between Air and Sea for a generic gas `X` with fixed Henry law coefficient `K_0`
and piston velocity `vpiston`.
The Reaction-local Variables `X_aq_conc`, `pXatm`, `X_airsea_exchange` should be linked to the appropriate variables using the
`variable_links:` section in the .yaml file.
# Parameters
$(PARS)
# Methods and Variables
$(METHODS_DO)
"""
Base.@kwdef mutable struct Reaction_Min_AirSeaExchange{P} <: PB.AbstractReaction

Expand All @@ -27,11 +37,10 @@ function PB.register_methods!(rj::Reaction_Min_AirSeaExchange)

vars = [

PB.VarDep( "CO2_aq_conc", "mol m-3", "CO2_aq concentration per cell"),
PB.VarDepScalar( "pCO2atm", "atm", "atmospheric pCO2, unit is atm"),
# PB.VarDepScalar( "pCO2PAL", "", "atmospheric pCO2 normalized to present day"),
PB.VarContrib( "CO2_airsea_exchange", "mol yr-1", "it is the calcalation for CO2_airsea_exchange"),
PB.VarDep( "area", "m2", "surface area"),
PB.VarDep( "X_aq_conc", "mol m-3", "ocean concentration per cell"),
PB.VarDepScalar( "pXatm", "atm", "atmospheric partial pressure, unit is atm"),
PB.VarContrib( "X_airsea_exchange", "mol yr-1", "calculated airsea exchange flux for gas X"),
PB.VarDep( "area", "m2", "surface area"),

]

Expand All @@ -43,18 +52,16 @@ end
# do method, called each main loop timestep
function do_Min_AirSeaExchange(m::PB.ReactionMethod, pars, (varsdata, ), cellrange::PB.AbstractCellRange, deltat)

rj = m.reaction

# mol m-3 mol m-3 atm-1 atm
CO2_aq_conc_eqb = pars.K_0[] * varsdata.pCO2atm[]
# mol m-3 mol m-3 atm-1 atm
X_aq_conc_eqb = pars.K_0[] * varsdata.pXatm[]

for i in cellrange.indices

# mol m-3 mol m-3
CO2_aq_conc = varsdata.CO2_aq_conc[i]
X_aq_conc = varsdata.X_aq_conc[i]

# mol yr-1 mol m-3 mol m-3 m yr-1 m2
varsdata.CO2_airsea_exchange[i] -= (CO2_aq_conc - CO2_aq_conc_eqb) * pars.vpiston[] * varsdata.area[i]
# mol yr-1 mol m-3 mol m-3 m yr-1 m2
varsdata.X_airsea_exchange[i] -= (X_aq_conc - X_aq_conc_eqb) * pars.vpiston[] * varsdata.area[i]

end

Expand Down
Loading

0 comments on commit 9268449

Please sign in to comment.