Skip to content

Commit

Permalink
Merge pull request #19 from PALEOtoolkit/ocean_chemistry
Browse files Browse the repository at this point in the history
Move Xiyuan's ocean chemistry example to a new folder
  • Loading branch information
sjdaines authored May 12, 2023
2 parents b2d1de8 + bd58c20 commit d1bc22f
Show file tree
Hide file tree
Showing 13 changed files with 494 additions and 177 deletions.
1 change: 1 addition & 0 deletions docs/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@ DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
PALEOboxes = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4"
PALEOmodel = "bf7b4fbe-ccb1-42c5-83c2-e6e9378b660c"
PALEOocean = "41de04b1-2efd-44ae-92ae-39d71a4fd99b"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
Expand Down
1 change: 1 addition & 0 deletions examples/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
PALEOboxes = "804b410e-d900-4b2a-9ecd-f5a06d4c1fd4"
PALEOcopse = "4a6ed817-0e58-48c6-8452-9e9afc8cb508"
PALEOmodel = "bf7b4fbe-ccb1-42c5-83c2-e6e9378b660c"
PALEOocean = "41de04b1-2efd-44ae-92ae-39d71a4fd99b"
PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
Expand Down
2 changes: 2 additions & 0 deletions examples/doc_order.txt
Original file line number Diff line number Diff line change
Expand Up @@ -2,3 +2,5 @@ reservoirs
solvers
error_examples
CPU
ocean_chemistry

136 changes: 136 additions & 0 deletions examples/ocean_chemistry/README.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
# Ocean chemistry

These examples illustrate how to implement:
- a minimal model of the marine carbonate system
- air-sea exchange of CO2 between ocean and atmosphere Domains

See `PALEOmodel` [documentation](https://paleotoolkit.github.io/PALEOmodel.jl/) for more information on analysing model output.

## 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).

### 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`:
```@eval
str = read("../../../../examples/ocean_chemistry/reactions_Alk_pH.jl", String)
str = """```julia
$str
```"""
import Markdown
Markdown.parse(str)
```

### 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.
```@eval
str = read("../../../../examples/ocean_chemistry/config_ex1.yaml", String)
str = """```julia
$str
```"""
import Markdown
Markdown.parse(str)
```

### Run script
The script to run the model (file `examples/ocean_chemistry/run_ex1.jl`) contains:
```@eval
str = read("../../../../examples/ocean_chemistry/run_ex1.jl", String)
str = """```julia
$str
```"""
import Markdown
Markdown.parse(str)
```
and produces output showing the change, if `TAlk_conc` increase, how the carbonic acid and pH change in ocean:
```@example ex1
include("../../../../examples/ocean_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
```

![](ex1_plot1.svg)
![](ex1_plot2.svg)
![](ex1_plot3.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 Air-sea exchange

Adds air-sea exchange of CO2 to Example 1

### 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.
```@eval
str = read("../../../../examples/ocean_chemistry/reactions_AirSeaExchange.jl", String)
str = """```julia
$str
```"""
import Markdown
Markdown.parse(str)
```

### 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
str = read("../../../../examples/ocean_chemistry/config_ex2.yaml", String)
str = """```julia
$str
```"""
import Markdown
Markdown.parse(str)
```

### Run script
The script to run the model (file `examples/ocean_chemistry/run_ex2.jl`) contains:
```@eval
str = read("../../../../examples/ocean_chemistry/run_ex2.jl", String)
str = """```julia
$str
```"""
import Markdown
Markdown.parse(str)
```
and produces output showing the change, if `TAlk_conc` increase, how the carbonic acid and pH change in ocean and CO2 change in the air:
```@example ex2
include("../../../../examples/ocean_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(plot(paleorun.output, "atm.pCO2atm", )) # hide
savefig("ex2_plot4.svg"); nothing # hide
display(plot(paleorun.output, ["global.C_total", "atm.CO2", "ocean.DIC_total"] ; ylabel="atm-ocean carbon (mol)"))
savefig("ex2_plot5.svg"); nothing # hide
```

![](ex2_plot1.svg)
![](ex2_plot2.svg)
![](ex2_plot3.svg)
![](ex2_plot4.svg)
![](ex2_plot5.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
```

For more information and cooperation, please communicate with us!
83 changes: 83 additions & 0 deletions examples/ocean_chemistry/config_ex1.yaml
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
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


solve_Alk_pH:

class: Reaction_Alk_pH

parameters:

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

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


oceanfloor:
# unused here, set up by ReactionOceanNoTransport

oceansurface:
# unused here, set up by ReactionOceanNoTransport
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
Minimal_Alk_pH:
Minimal_Alk_pH_AirSea:

domains:

Expand All @@ -7,18 +7,31 @@ Minimal_Alk_pH:

reactions:

sum_E:
sum_C:

class: ReactionSum

parameters:

vars_to_add: ["ocean.DIC", "atm.CO2"]
vars_to_add: ["ocean.DIC_total", "atm.CO2"]

variable_links:

sum: C_total

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 !!



# --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Expand All @@ -40,7 +53,7 @@ Minimal_Alk_pH:

reservoir_TA:

class: ReactionReservoir
class: ReactionReservoirTotal

variable_links:
R*: TAlk*
Expand All @@ -51,7 +64,7 @@ Minimal_Alk_pH:


reservoir_DIC:
class: ReactionReservoir
class: ReactionReservoirTotal

variable_links:
R*: DIC*
Expand All @@ -60,22 +73,25 @@ Minimal_Alk_pH:
R:norm_value: 2.0e0



solve_Alk_pH:

class: Reaction_Alk_pH

parameters:

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

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
B_total: 4.2e-4

variable_links:
density: rho_ref # OceanNoTransport provides density as rho_ref


# --------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Expand All @@ -84,9 +100,27 @@ Minimal_Alk_pH:
# air-sea exchange reactions

reactions:
# calculate air-sea exchange flux
solve_AirSea_Exchange:

class: Reaction_Min_AirSeaExchange

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

variable_links:

pCO2atm: atm.pCO2atm
CO2_aq_conc: ocean.oceansurface.CO2_aq_conc

area: Asurf # ocean surface area

CO2_airsea_exchange: fluxAtmtoOceansurface.flux_CO2

# apply air-sea flux to ocean surface

# apply air-sea flux to ocean surface
transfer_fluxAtmtoOceansurface:

class: ReactionFluxTransfer
Expand Down Expand Up @@ -129,31 +163,9 @@ Minimal_Alk_pH:

variable_attributes:

R:norm_value: 4.956e16
R:norm_value: 4.956e16 # 280e-6 ppm
R:initial_value: 4.956e16


solve_AirSea_Exchange:

class: Reaction_Min_AirSeaExchange

parameters:

K_0: 3.4e-5
vpiston: 1138.8

variable_links:


CO2_aq_conc: ocean.CO2_aq_conc

area: ocean.Abox

CO2_airsea_exchange: fluxAtmtoOceansurface.flux_CO2




transfer_AtmtoOceansurface:

class: ReactionFluxTransfer
Expand Down
Loading

0 comments on commit d1bc22f

Please sign in to comment.