diff --git a/docs/src/paleo_references.bib b/docs/src/paleo_references.bib index 47eabc1..835041b 100644 --- a/docs/src/paleo_references.bib +++ b/docs/src/paleo_references.bib @@ -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}, diff --git a/examples/doc_order.txt b/examples/doc_order.txt index 5c7d363..cdf1894 100644 --- a/examples/doc_order.txt +++ b/examples/doc_order.txt @@ -2,5 +2,6 @@ reservoirs solvers error_examples CPU +CPU_modular ocean_chemistry diff --git a/examples/ocean_chemistry/README.md b/examples/ocean_chemistry/README.md index 248f1a3..a6c9988 100644 --- a/examples/ocean_chemistry/README.md +++ b/examples/ocean_chemistry/README.md @@ -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 @@ -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. @@ -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 @@ -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 @@ -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 diff --git a/examples/ocean_chemistry/config_ex1.yaml b/examples/ocean_chemistry/config_ex1.yaml index dab50f4..3faac3c 100644 --- a/examples/ocean_chemistry/config_ex1.yaml +++ b/examples/ocean_chemistry/config_ex1.yaml @@ -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 diff --git a/examples/ocean_chemistry/config_ex2.yaml b/examples/ocean_chemistry/config_ex2.yaml index 3f561da..00f7d25 100644 --- a/examples/ocean_chemistry/config_ex2.yaml +++ b/examples/ocean_chemistry/config_ex2.yaml @@ -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 @@ -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 diff --git a/examples/ocean_chemistry/reactions_AirSeaExchange.jl b/examples/ocean_chemistry/reactions_AirSeaExchange.jl index 1479fe7..350b4f5 100644 --- a/examples/ocean_chemistry/reactions_AirSeaExchange.jl +++ b/examples/ocean_chemistry/reactions_AirSeaExchange.jl @@ -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 @@ -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"), ] @@ -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 diff --git a/examples/ocean_chemistry/reactions_Alk_pH.jl b/examples/ocean_chemistry/reactions_Alk_pH.jl index d72e0e2..3cf8b10 100644 --- a/examples/ocean_chemistry/reactions_Alk_pH.jl +++ b/examples/ocean_chemistry/reactions_Alk_pH.jl @@ -1,51 +1,52 @@ module Min_Alk_pH import PALEOboxes as PB +using PALEOboxes.DocStrings # for $(PARS) and $(METHODS_DO) """ - Min_Alk_pH + Reaction_Alk_pH -Minimal example, TAlk_conc is first order decay of a variable. +Minimal example for aqueous carbonate system. +Solves for carbon, boron and water species given `pH`, calculates difference `TAlk_error` from +required alkalinity. +Use in conjunction with a DAE solver, where this Reaction provides an algebraic constraint `TAlk_error` on `pH`. + +# Parameters +$(PARS) + +# Methods and Variables +$(METHODS_DO) """ Base.@kwdef mutable struct Reaction_Alk_pH{P} <: PB.AbstractReaction base::PB.ReactionBase pars::P = PB.ParametersTuple( - # PB.ParDouble("kappa", 1.0, units="yr-1", description="first order decay constant"), - # PB.ParDouble("TAlk_conc_change", 10.0e-6, units="mol kg-1", description="TA concentration change every kappa"), - # PB.ParDouble("density", 1027.0, units="kg m-3", description="current seawater density"), PB.ParDouble("K_1", 1.4e-6, units="mol kg-1", description="equilibrium constant of CO2_aq and HCO3-"), PB.ParDouble("K_2", 1.2e-9, units="mol kg-1", description="equilibrium constant of HCO3- and CO32-"), - PB.ParDouble("K_w", 6.0e-14, units="mol2 kg-2", description="equilibrium constant of water at S=35, T=25°C"), - PB.ParDouble("K_B", 2.5e-9, units="mol kg-1", description="equilibrium constant of B(OH)4-"), - PB.ParDouble("B_total", 4.2e-4, units="mol kg-1", description="total concentrations of B(OH)4- and B(OH)3"), - + PB.ParDouble("K_w", 6.0e-14, units="mol^2 kg-2", description="equilibrium constant of water at S=35, T=25°C"), + PB.ParDouble("K_B", 2.5e-9, units="mol kg-1", description="equilibrium constant of B(OH)4-"), ) end function PB.register_methods!(rj::Reaction_Alk_pH) vars = [ - # PB.VarDep("DIC", "mol", "reservoir for species DIC"), PB.VarDep("DIC_conc", "mol m-3", "DIC concentration"), - # PB.VarContrib("DIC_sms", "mol yr-1", "reservoir DIC source - sink"), - # PB.VarDep("TAlk", "mol", "reservoir for species TA"), - PB.VarDep("TAlk_conc", "mol m-3", "TA concentration"), - # PB.VarContrib("TAlk_sms", "mol yr-1", "reservoir TA source - sink"), - # PB.VarProp("TAlk_decay_flux", "mol yr-1", "decay flux from reservoir TA"), + PB.VarDep("TAlk_conc", "mol m-3", "TA concentration"), + PB.VarDep("B_conc", "mol m-3", "total Boron concentration"), PB.VarConstraint("TAlk_error", "mol m-3", "in order to solve TA, we set it"), PB.VarProp("HCO3_conc", "mol m-3", "HCO3- concentration"), PB.VarProp("CO3_conc", "mol m-3", "CO32- concentration"), PB.VarProp("CO2_aq_conc", "mol m-3", "CO2_aq concentration"), + PB.VarProp("BOH3_conc", "mol m-3", "BOH3 concentration"), PB.VarProp("BOH4_conc", "mol m-3", "BOH4- concentration"), PB.VarProp("H_conc", "mol m-3", "concentration of H+"), PB.VarProp("OH_conc", "mol m-3", "concentration of OH-"), PB.VarState("pH", "", "it is the calcalation for pH"), - # PB.VarDep("volume", "m3", "ocean volume"), - PB.VarDep("density", "kg m-3", "ocean density"), + PB.VarDep("density", "kg m-3", "ocean density"), ] PB.add_method_do!(rj, do_Min_Alk_pH, (PB.VarList_namedtuple(vars), ) ) @@ -55,6 +56,8 @@ function PB.register_methods!(rj::Reaction_Alk_pH) return nothing end +# called at model start +# provide an initial value for pH (exact value is not important, just needs to be in a reasonable range) function setup_carbchem( m::PB.ReactionMethod, pars, (vars, ), cellrange::PB.AbstractCellRange, attribute_name ) attribute_name in (:initial_value, :norm_value) || return @@ -79,22 +82,10 @@ function do_Min_Alk_pH(m::PB.ReactionMethod, pars, (varsdata, ), cellrange::PB.A for i in cellrange.indices density = varsdata.density[i] - # mol yr-1 yr-1 mol kg-1 kg m-3 m3 - # varsdata.TAlk_decay_flux[i] = pars.kappa[] * pars.TAlk_conc_change[] * density * varsdata.volume[i] - - # mol yr-1 mol yr-1 - # varsdata.TAlk_sms[i] = varsdata.TAlk_decay_flux[i] - # mol kg-1 mol m-3 kg m-3 - DIC_conc_kg = varsdata.DIC_conc[i] / density - - - # mol kg-1 mol m-3 kg m-3 - # TAlk_conc_kg = varsdata.TAlk_conc[i] / density + DIC_conc_kg = varsdata.DIC_conc[i] / density + B_total_kg = varsdata.B_conc[i] / density - # mol kg-1 = mol kg-1 - # H_conc_kg = 10^( -varsdata.pH[i] ) - # SD units ... # mol kg-1 = mol L-1 / kg m-3 * L m-3 H_conc_kg = 10^( -varsdata.pH[i] ) / density * 1000.0 @@ -106,9 +97,11 @@ function do_Min_Alk_pH(m::PB.ReactionMethod, pars, (varsdata, ), cellrange::PB.A # mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1 mol kg-1 CO3_conc_kg = DIC_conc_kg / ( 1 + H_conc_kg/pars.K_2[] + ((H_conc_kg)^2)/(pars.K_1[]*pars.K_2[]) ) - + # mol kg-1 mol kg-1 mol kg-1 mol kg-1 - BOH4_conc_kg = pars.B_total[] / ( 1 + H_conc_kg/pars.K_B[] ) + BOH4_conc_kg = B_total_kg / ( 1 + H_conc_kg/pars.K_B[] ) + + BOH3_conc_kg = B_total_kg - BOH4_conc_kg # mol kg-1 mol2 kg-2 mol kg-1 OH_conc_kg = pars.K_w[] / H_conc_kg @@ -118,21 +111,12 @@ function do_Min_Alk_pH(m::PB.ReactionMethod, pars, (varsdata, ), cellrange::PB.A # mol m-3 kg m-3 mol kg-1 varsdata.H_conc[i] = density * H_conc_kg - - # mol m-3 kg m-3 mol kg-1 varsdata.OH_conc[i] = density * OH_conc_kg - - # mol m-3 kg m-3 mol kg-1 varsdata.HCO3_conc[i] = density * HCO3_conc_kg - - # mol m-3 kg m-3 mol kg-1 varsdata.CO3_conc[i] = density * CO3_conc_kg - - # mol m-3 kg m-3 mol kg-1 varsdata.CO2_aq_conc[i] = density * CO2_aq_conc_kg - - # mol m-3 kg m-3 mol kg-1 - varsdata.BOH4_conc[i] = density * BOH4_conc_kg + varsdata.BOH4_conc[i] = density * BOH4_conc_kg + varsdata.BOH3_conc[i] = density * BOH3_conc_kg # mol m-3 mol m-3 mol kg-1 kg m-3 varsdata.TAlk_error[i] = varsdata.TAlk_conc[i] - TAlk_conc_kg_calcu * density diff --git a/examples/ocean_chemistry/run_ex1.jl b/examples/ocean_chemistry/run_ex1.jl index 2bd6d47..fb13e17 100644 --- a/examples/ocean_chemistry/run_ex1.jl +++ b/examples/ocean_chemistry/run_ex1.jl @@ -51,6 +51,21 @@ display(plot(paleorun.output, ["ocean.TAlk_conc", "ocean.DIC_conc"], 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)")) -# display(plot(paleorun.output, "ocean.DIC", xlims=(0, 150.0), (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, + ) +) +